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SUMMARY 

Using  invariance  principles  in  a  systematic  fashion,  we 
shall  derive  not  only  new  analytic  formulations  of  the 
classical  partiole  processes,  those  of  transport  theory, 
radiative  transfer;  random  walk,  multiple  scattering,  and 
diffusion  theory,  but,  in  addition,  new  computational 
algorithms  which  seem  well  fitted  to  the  capabilities  of 
digital  computers.  Whereas  the  usual  methods  reduce  problems 
to  the  solution  of  systems  of  linear  equations,  we  shall  try 
to  reduce  problems  to  the  Iteration  of  nonlinear  transformations . 

Although  we  have  analogous  formulations  of  wave  processes, 
we  shall  reserve  for  a  second  paper  in  this  series  a  detailed 
and  extensive  treatment  of  this  part  of  mathematical  physics. 
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INVARIANT  IMBEDDING  AND  MATHEMATICAL  PHYSICS— I: 

PARTICLE  PROCESSES 

Richard  Bellman,  Robert  Kalaba  and  G.  Milton  Wing 

Part  I 

INTRODUCTION 

The  classical  equations  of  mathematical  physics  can  be 
put  In  the  form 

(1)  ut  -  T(u) 

where  u  Is  a  vector  function  of  a  space  vector  p,  restricted 
to  a  region  R,  and  the  time,  u  =»  u(p,t).  The  operator  T  la 
In  many  cases  a  linear  partial  differential  operator,  In  some 
cases  a  linear  Integral  operator,  and.  If  we  Insist  upon 
realism,  a  nonlinear  operator.  The  steady  state  version  is 
obtained  by  setting  the  vector  ut  equal  to  zero. 

Since  equations  of  this  type  usually  have  an  infinite 
number  of  solutions,  it  is  necessary  to  attach  some  further 
restrictions  in  order  to  single  out  a  particular  solution.  To 
do  this,  we  usually  assign  initial  values, 

(2)  u(p,0)  =  v ( p ) ,  p  c  R, 
and  boundary  valuee 

(3)  u(p,t)  =  w(p,t), 

for  p  €  B,  the  boundary  of  R. 
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Problems  of  this  nature  have  two  types  of  difficulties 
associated,  difficulties  which  are,  of  course,  inseparably 
intertwined,  those  of  analytic  character  and  those  of  compu¬ 
tational  nature.  Among  the  many  methods  which  have  been  pro¬ 
posed  is  the  theory  of  semigroups.  The  guiding  ideas  were 
first  enunciated  by  Hadam&rd,  and  subsequently  were 
systematically  pursued  oy  HI lie  and  Yoslda;  see  HI  lie  and 
Phillips,  M.  for  a  thorough  exposition  and  many  references. 
Classically,  the  semigroup  concept  has  been  exploited  in  the 
time  domain.  Our  aim  is  to  show  that  this  basic  method  can  be 
applied  in  a  much  wider  area,  using  other  physical  variables 
of  significance  as  semigroup  variables. 

Using  l-nvariance  principles  in  a  systematic  fashion,  we 
shall  derive  not  only  new  analytic  formulations  of  the  classi¬ 
cal  particle  processes,  those  of  transport  theory,  radiative 
transfer,  random  walk,  multiple  scattering,  and  diffusion 
theory,  but,  in  addition,  new  computational  algorithms  which 
seem  well  fitted  to  the  capabilities  of  digital  computers. 
Whereas  the  usual  methods  reduce  problems  to  the  solution  of 
systems  of  linear  equations,  we  shall  try  to  reduce  problems 
N  the  Iteration  of  nonlinear  transformations. 

Although  we  have  analogous  formulations  of  wave  processes, 
briefly  presented  in  [2, 3, 4],  we  shall  reserve  for  a  second 
paper  in  this  series  a  detailed  and  extensive  treatment  of 
this  part  of  mathematical  physics. 

Our  interest  in  the  field  of  invariance  principles  was 
aroused  by  the  elegant  and  fundamental  work  of  Chandrasekhar 
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In  the  theory  of  radiative  transfer  [30]  .  His  results,  in  turn, 
are  generalizations  of  those  of  Ambarzuialan  who  seems  to  have 
been  the  first  to  have  consciously  employed  invariance  principles 
in  any  significant  fashion  £29].  Since  then,  in  addition 
to  our  work,  reference  to  which  nay  be  found  at  the  end  of  the 
paper,  there  have  been  important  contributions  by  Preiaendorfer, 
[5,6],  Ramakrishnan,  [7],  Redheffer,  [8],  and  Ueno,  Q?,10,ll]  . 

In  addition,  there  are  some  unpublished  results  due  to  Harris. 

Independently,  functional  equation  techniques  were  intro¬ 
duced  into  the  theory  of  branching  processes,  in  particular, 
those  arising  in  cosmic  ray  cascade  theory  and  biological 
mutation,  by  Bellman  and  Harris,  [12,15] ,  and  Janossy,  [l4]  . 
Surveys  of  the  many  results  obtained  over  the  last  ten  years 
may  be  found  in  Harris,  [15,16],  and  Ramakrishnan,  [17]  .  Also 
in  the  theory  of  dynamic  programming,  [18]  ,  in  connection  with 
the  treatment  of  minimization  and  maximization  problems,  we 
find  a  use  of  invariance  principles  and  functional  equations 
which  is  quite  similar  in  spfc'lt  to  what  we  shall  find  below  in 
the  treatment  of  purely  descriptive  processes. 

In  plaoe  of  beginning  with  an  abstract  formulation  of 
particle  processes  and  an  abstract  presentation  of  the 
principles  of  invariant  imbedding,  we  shall  start  with  a  study 
of  a  particular  process,  neutron  transport  and  multiplication. 

The  difference  in  formulation  between  the  usual  approach  and 
that  furnished  by  "invariant  imbedding,"  as  we  shall  call  our 
systematic  application  of  invariance  principles,  will  readily 
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be  seen.  Nevertheless,  as  we  shall  show,  both  are  merely 
particular  Instances  of  a  general  approach. 

Having  gone  through  a  spectrum  of  transport  processes, 
steady-state  and  time-dependent,  energy- Independent  and 
energy-dependent,  one-dlmenslonal  and  multidimensional, 
unchanging  medium  and  Stefan- type,  we  shall  abstract  the 
basic  Ideas  of  invariant  Imbedding. 

Following  this,  we  shall  apply  these  techniques  to  the 
study  of  random  walk  and  multiple  scattering,  to  the  study  of 
radiative  transfer  and  diffusion.  Our  treatment  of  these 
fields  will  be  much  briefer  since  much  of  what  is  done  in  the 
part  devoted  to  neutron  transport  can  easily  be  transcribed 
and  applied  in  these  other  areas. 

In  what  follows,  we  shall  pursue  a  purely  formal  path, 
leaving  aside  all  questions  of  existence,  uniqueness  and  so 
on.  What  is  interesting,  however,  is  that  our  approach 
enables  us  to  handle  many  of  these  questions  in  a  much  simpler 
and  straightforward  way  than  that  furnished  by  the  conventional 
road . 

Although  we  are  in  part  motivated  by  a  search  for 
feasible  computational  techniques,  we  shall  actually  avoid  any 
discussion  of  actual  numerical  techniques.  In  subsequent 
papers,  we  shall  treat  these  matters  in  great  detail.  Here  we 
shall  restrain  ourselves  to  generalities. 

The  equations  of  invariant  imbedding  are  related  to  the 
variational  formulas  of  Hadamard  type,  expressing  the 
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dependence  of  the  Green's  function  of  a  region  upon  the 
dimensions  of  the  region  (see  [52]). 

Finally,  let  us  note  that  no  previous  knowledge  of  the 
equations  of  mathematical  physics  is  required.  All  equations 
will  be  derived  from  first  principles,  directly  from  the 
mathematical  model  of  the  physical  process. 
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Part  II 

NEUTRON  TRANSPORT  AND  MULTIPLICATION 

1.  Introduction 

Let  us  begin  our  journey  with  the  examination  of  a  num¬ 
ber  of  intriguing  mathematical  problems  which  arise  in  the 
study  of  various  aspects  of  neutron  transport  and  multiplica¬ 
tion.  A  consideration  of  some  of  the  many  different  hypotheses 
that  can  be  made  will  give  us  an  opportunity  to  display  the 
versatility  of  the  theory  of  invariant  imbedding. 

Our  basic  assumption  is  that  a  neutron  is  a  point  particle 
which  is  completely  specified  at  any  time  by  its  direction  of 
motion  and  its  energy.  These  two  properties  determine  its 
state.  As  the  neutron  traverses  the  medium  within  which  the 
transport  process  takes  place,  it  suffers  certain  changes  of 
state,  i.e.,  changes  in  energy  and  direction,  due  to  inter¬ 
actions  with  the  medium  and  with  other  neutrons.  In  addition, 
we  have  the  relatively  new  and  very  Important  phenomena  of 
fission.  Certain  interactions  can  result  in  an  increase  in 
the  number  of  neutrons,  the  fission  process. 

The  probabilities  of  these  events  are  measured  by  "cross- 
sections”  or  "mean  free  paths."  Occasionally,  we  shall  talk 
in  deterministic  terms,  and  occasionally  in  stochastic  terms, 
depending  upon  which  is  more  convenient.  The  difference  is 
mor*  apparent  than  real,  since  the  use  of  expected  values  in  a 
stochastic  model  leads  to  a  completely  deterministic  version 
based  upon  fluxes. 
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Within  the  framework  of  a  mathematical  model  conetructed 
along  these  lines,  a  model  we  shall  make  more  precise  In  the 
following  section,  we  wish  to  explain  and  predict  the  pheno¬ 
menon  of  criticality,  and  to  determine  the  internal  and 
external  fluxes  as  functions  of  the  spatial  dimensions,  the 
time  and  other  parameters.  Problems  of  this  nature  are  of 
great  complexity  from  the  mathematical  side,  and  thus  of  even 
greater  fascination,  even  when  greatly  simplified  physical 
models  are  used.  When  more  realistic  assumptions  are  made, 
the  analytic  aspects  become  truly  formidable,  and  the  problem 
of  obtaining  numerical  results  much  more  burdensome. 

It  is  not  to  be  expected  that  recondite  scientific 
questions  will  yield  readily  to  any  single  approach.  Rather 
it  is  to  be  expected  that  with  the  aid  of  a  battery  of  methods, 
each  of  which  chips  away  at  some  of  the  obstacles  in  our  path, 
we  can  eventually  clear  a  road  which  will  take  us  some  distance 
toward  our  goal. 

The  classical  equations  of  transport  theory  can  be 

1 

effectively  applied  in  a  number  of  cases.  Approximate  methods 
of  various  degrees  of  efficacy  and  associated  results  may  be 
found  in  the  book  by  Davison  [19] .  Rigorous  discussion  of  these 
techniques  can  lead  to  quite  complex  analysis;  see  for  example 
the  papers  by  Lehner  and  Wing  [20,21,44,49] ,  Jorgens  [43],  and 
Pimbley  [48 j . 

A  number  of  questions  can  be  Btudied  by  means  of  the 
mathematical  theory  of  branching  processes.  The  study  of  age- 
independent  processes  was  begun  by  Harris  [ 15,16],  and 
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Everett  and  Ulan,  [22] ,  Independently  of  each  other. 

Eeeentlally  it  reduces  to  the  study  of  the  iteration  of  power 
series,  with  probabilistic  overtones.  The  theory  of  age- 
dependent  branching  processes,  based  upon  the  systematic  usage 
of  functional  equations,  was  begun  by  Bellman  and  Harris, 

[12,13],  and  independently  by  Janossy,  [14],  Detailed  expositions 
with  many  references  will  be  found  in  the  monograph  by  Harris, 
[l^J,  and  the  expository  papers  by  Harris,  [l^,  and 
Ramakrishnan,  [17]  . 

It  is  natural  to  construct  simplified  models  in  a  situation 
characterized  by  severe  mathematical  difficulties  and  by 
physical  complexity  as  well.  The  usual  hope  is  that  the 
exploration  of  these  models  will  furnish  valuable  experience 
and  that  the  understanding  of  these  more  transparent  models 
will  enable  us  to  penetrate  into  the  more  obscure  versions. 
However,  as  mentioned  above,  even  apparently  simple  processes 
give  rise  to  sophisticated  analysis. 

Furthermore,  as  we  shall  discuss  repeatedly  below,  unless 
the  problems  are  carefully  formulated  they  cannot  be  resolved 
in  numerical  terms  in  any  straightforward  fashion.  Our  objec¬ 
tive  in  the  pages  that  follow  is  to  formulate  a  variety  of 
transport  processes  in  a  way  which  will  permit  us  to  obtain 
numerical  solutions  with  the  aid  of  digital  computers.  As  is 
often  the  case  In  mathematics  and  physics,  a  significant 
Improvement  in  computational  technique  requires  a  new  conceptual 
and  analytic  approach. 
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It  turns  out  that  In  the  process  of  fulfilling  one  of 
our  goals,  numerical  solution  of  problems,  we  obtain 
as  byproducts  a  host  of  Interesting  and  elegant  analytic  re¬ 
sults,  together  with  powerful  methods  for  establishing 
existence  and  uniqueness  theorems  for  the  associated  functional 
equations  and  for  the  classical  functional  equations  of  mathe¬ 
matical  physics.  Many  of  these  equations  are  quite  difficult 
to  treat  along  conventional  lines. 

Our  principal  tool  will  be  the  theory  of  invariant  imbed¬ 
ding.  Rather  than  attempt  to  define  precisely  what  turns  out 
to  be  more  a  state  of  mind  than  anything  else,  we  shall  first 
give  a  number  of  applications  of  the  methods.  Subsequently, 
we  shall  try  to  distill  the  essence  of  these. 

2.  A  Simple  Neutron  Transport  and  Multiplication  Process 

Let  us  now  describe  a  simple  mathematical  model  of  a 
neutron  transport  process  with  fission.  Subsequently,  we 
shall  add  a  number  of  interesting  features  such  as  collision 
between  neutrons,  energy  and  time  dependence,  and  so  on.  For 
the  immediate  purpose  of  illustrating  both  how  the  classical 
approach  is  made  and  how  invariant  imbedding  techniques  are 
used,  there  are  great  advantages  to  using  the  simplest  possible 
version  possessing  certain  structural  properties. 

As  noted  above,  we  take  the  neutron  to  be  a  point  particle, 
and  we  allow  at  the  moment  only  one- dimensional  motion  along  a 
line,  or  part  of  a  line.  To  simplify  matters  still  further,  we 
assume  that  there  is  no  energy  dependence.  As  this  blithe. 
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carefree  neutron  inodes  along  the  line,  it  may  suffer  a  collision 
with  the  constituent  elements  of  the  line.  Again  to  simplify 
the  algebra,  we  suppose  that  only  fission  collisions  occur, 
resulting  in  one  neutron  moving  to  the  left  and  one  to  the 
right.  This  is  the  only  type  of  interaction  we  shall  allow 
between  the  neutron  and  the  transport  medium  at  the  moment. 
Furthermore,  we  shall  suppose  that  there  are  no  neutron-neutron 
interactions . 

To  make  this  verbiage  precise,  let  us  consider  a  finite 
interval  [0,x]  (the  reason  for  this  apparently  loose  usage 
of  x  to  designate  an  endpoint  will  be  made  clear  subsequently. 
At  the  present,  let  us  merely  state  that  it  is  done  with  malice 
aforethought.),  a  one-dimensional  rod,  with  the  following 
properties: 

(l)  (a)  When  a  neutron  traverses  an  infinitesimal  length  A, 

in  either  direction,  there  is  a  probability  o A  +  o(a) 
that  fission  will  occur. 

(b)  When  fission  occurs,  two  neutrons  are  produced,  one 
going  to  the  right  and  one  to  the  left-  Each  of 
these  has  the  same  properties  as  the  original  neutron. 


0  y  y+A  x 

Fig.  1 

The  notation  f(x)  -  o(g(x))  is  used  to  mean 
lim  f(x)/g(x)  ■  0,  where  the  sense  of  the  limit  is  usually 
obvious  from  the  context. 
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(c)  There  Is  a  probability  1  -  oA  +  o(A)  that  no  inter¬ 
action  occurs  in  A ,  which  means  no  change  in  the 
direction  of  the  neutron. 

(d)  When  a  neutron  leaves  the  rod,  it  cannot  return  and  it 
has  no  further  effect  upon  the  transport  process. 

It  would  not  be  difficult  to  include  absorption  effects 
and  collisions  which  merely  change  direction  or  to  allow  R 
neutrons,  R  /  2,  out  of  a  collision.  Since  these  effects  are 
treated  subsequently,  in  multiple  scattering,  radiative  transfer 
and  random  walk,  we  shall  omit  them  here  to  keep  the  analytic 
details  to  an  irreducible  minimum. 

The  quantity  o  is  called  the  "macroscopic  cross-section." 
Occasionally,  we  shall  write  it  as  l/A,  where  7\  is  called 
the  "mean  free  path."  If  the  rod  is  homogeneous,  these  quanti¬ 
ties  sre  constant,  otherwise  we  write  c(y)  and  7\{y)  for 
the  quantities  associated  with  the  interval  [y,y  +  a]. 

We  shall  begin  by  considering  steady-state  neutron  flux. 

The  more  general  time-dependent  case  will  be  considered  below. 
Let  a  unit  flux  of  neutrons  (that  is,  one  neutron  per  unit  time) 
be  incident  upon  the  right  end  of  the  rod,  and  let  it  be  desired 
to  determine  the  right  and  left  fluxes  at  any  internal  point  y, 
as  well  as  the  fluxes  out  at  zero  and  x.  These  we  shall  regul¬ 
arly  refer  to  as  transmitted  and  reflected  fluxes,  respectively. 

Our  first  formulation  will  be  the  classical  one,  resulting 
in  simple  versions  of  the  linearized  Boltzmann  equation. 
Introduce  the  functions 
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(2)  uR(y)  °  the  expected  number  of  neutrons  going  to  the 

right  at  y  per  unit  time, 

uL(y)  -  the  expected  number  of  neutrons  going  to  the 
left  at  y  per  unit  time. 

To  obtain  differential  equations  for  uR  and  u^>  we 
apply  simple  conservation  laws  for  the  right  and  left-hand 
flows  at  y.  These  are  input-output  equations  expressing  the 
fact  that  what  goes  out  is  the  sum  of  what  comes  in  and  what 
is  produced.  By  virtue  of  our  assumptions  concerning  the 
transport  and  fission  process  and  the  elementary  laws  of 
probability,  we  obtain  the  equations 

(3)  uR(y)  -  uR(y  -  A ) ( 1  -  oA)  +  (uR(y)  4-  uL(y))oA  +  o(a), 

uL(y)  =  uL(y  +  A )( 1  -  aA)  +  (uR(y)  +  uL(y))aA  +  o(a). 

Passing  to  the  limit  as  A  — *  0,  we  obtain  the  system  or 
differential  equations 

(M  uR ' ( y )  -  ouL(y), 

•  uL* (y)  -  ~  ouR(y). 

The  boundary  conditions  are 

uL(x)  =  1, 

ur(0)  -  0. 


(5) 
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These  express  the  fact  that  there  is  an  incident  flux  of 
unit  strength  at  the  point  x,  and  the  fact  that  there  is  no 
incident  flux  at  the  point  0.  Observe  a  property  which  we 
shall  repeatedly  stress:  The  physical  process  automatically 
leads  to  a  two-point  boundary  value  problem  when  formulated 
In  the  foregoing  way.  The  reason  for  discussing  this  fact 
will  be  discussed  in  detail  below. 

Finally,  let  us  note  that  we  can  obtain  the  most  general 
second  order  Sturm- Liouville  equation  from  the  foregoing 
process  i f  we  assume  that  right-hand  motion  at  y  has  a 
different  mean- free  path  than  left-hand  motion  and  that 
these  mean  free  paths  vary  with  y. 

3.  Invariant  Imbedding  Approach— Metaphysical 

We  now  wish  to  formulate  the  transport  process  described 
in  the  foregoing  section  in  different  terms.  Our  approach 
will  be  based  upon  the  theory  of  Invariant  imbedding.  What 
we  wish  to  do  is  to  imbed  the  particular  process  considered 
above  within  a  family  of  processes  of  similar  nature.  Al¬ 
though  this  appears  to  complicate  rather  than  simplify  the 
problem,  its  Justification  lies  in  the  fact  that  there  will 
exist  simple  relations  between  various  members  of  the  family 
which  can  be  utilised  to  determine  the  characteristics  of  a 
particular  member  cf  the  family. 

The  fact  that  the  structure,  or  anatomy,  of  a  particular 
organism  can  be  understood  quite  readily  in  terms  of  the  com- 
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paratlve  anatomy  of  a  phylum  la  well  established  in  the  field 
of  biology.  In  chemistry,  the  construction  of  the  Hendelieff- 
Moseley  periodic  table  of  the  elements  was  a  decisive  step 
forward.  In  mathematics,  the  method  of  continuity  is  one  of 
the  basic  devices  of  analysis  and  geometry.  It  follows  that 
in  pursuing  this  approach,  we  are  invoking  a  factotum  of 
science. 

Consider  the  way  in  which  an  experimental  physicist  might 
study  neutron  flux  in  a  rod.  Starting  with  a  rod  of  fixed 
length,  he  would  measure  reflected  and  transmitted  fluxes. 
Increasing  or  decreasing  the  length,  measurement  would  be  made 
of  the  corresponding  quantities.  The  final  data  would  consist 
of  two  curves,  one  the  reflected  flux  as  a  function  of  the 
length  of  the  rod,  the  other  the  transmitted  flux.  These 
would  be  functions  of  x,  the  length  of  the  rod. 

Our  aim  here  is  to  carry  out  the  analytic  equivalent  of 
this  program.  In  order  for  these  concepts  to  be  meaningful, 
we  must  find  a  way  of  relating  the  reflected  and  transmitted 
flux  for  a  rod  of  length  x  with  the  corresponding  fluxes  for 
rods  of  different  length.  We  propose  then  to  consider  the  set 
of  processes  obtained  by  letting  x  assume  any  positive 
value.  Our  choice  of  the  symbol  x  obviously  presages  this 
development. 

One  advantage  of  this  approach  as  far  as  reduction  of 
data  is  concerned  is  that  it  permits  a  direct  comparison  of 
analytic  results  with  experimental  results.  The  analytic  and 
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computational  advantages  will  be  discussed  in  extenso  below 
after  we  have  supplied  some  analytic  content  to  this  meta¬ 
physical  discourse. 


4.  Invariant  Imbedding  Approach — Analytical 
We  begin  by  introducing  the  function 


(l)  u(x)  -  the  expected  number  of  neutrons  reflected 

from  [0,x]  per  unit  time  as  a  result  of  an 
incident  unit  flux  of  neutrons  per  unit  time 
at  x. 


Transmitted  flux 


0 


- *  Reflected 

Plux 

—| -  4 -  Incident 

x-A  x  Plux 


Fig.  2 


Let  us  take  A  to  be  an  infinitesimal.  As  the  incident 
flux  passes  through  the  segment  ^x  -  A,x],  some  of  the 
neutrons  cause  fission  and  others  pass  through  unaffected  to 
become  incident  upon  [o,x  —  a|.  When  a  fission  occurs  in  A 
one  fission  neutron  emerges  at  x,  while  the  other  becomes  a 
part  of  the  Incident  flux  at  x  —  A. 

Some  of  the  neutrons  reflected  from  j^0,x  -  a]  may  cause 
flS3ion  while  passing  through  ^x  -  A,xj.  The  products  of  this 
fission  yield  a  contribution  to  the  reflected  flux  at  x  and 
furnish  another  source  of  neutrons  Incident  upon  [o,x  —  a|  . 

Fortunately,  although  the  physical  process  and  mathematical 
counterpart  are  exceedingly  complex  if  account  of  all  fissions 
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and  reflections  is  taken,  this  intricate  bookkeeping  is 

unnecessary  if  A  is  an  infinitesimal.  All  other  events, 

apart  from  those  taken  account  of  above,  have  a  probability 

2 

of  occurrence  of  order  A  or  higher.  Hence,  they  can  be 
neglected  in  the  derivation  of  the  differential  equation  for 
the  expected  flux  u(x).  Adding  up  the  various  effects  and 
their  associated  probabilities,  we  obtain  the  equation 

(2)  u(x)  -  oA(l  +  u(x  -  A)) 

+  (1  -  oA)[u(x  —  A)((l  -  oa)  +  oa(1  +  u(x  -  A))]]  +  o(A). 
Letting  A  — ♦  0,  we  derive  the  differential  equation 

(3)  u'(x)  •  d(l  ♦  u2(x)),  u(0)  -  0. 

This  first  order  nonlinear  differential  equation  is  called  a 
Riccati  equation.  As  we  shall  see,  this  type  of  quadratically 
nonlinear  equation  is  characteristic  of  the  equations  derived 
by  invariant  Imbedding  techniques.  In  contrast,  the  classical 
equations  are  linear.  Since  we  are  describing  the  same  pro¬ 
cess  in  different  ways,  there  must  be  relations  between  the 
analytic  descriptions.  We  shall  obtain  these  below. 

A  further  useful  function  is 

(4)  v(x)  «■  the  expected  transmitted  flux  per  unit  time 

as  a  result  of  a  unit  flux  per  unit  time 
incident  at  x. 

The  same  reasoning  as  above  shows  that  v(x)  satisfies  the 
equation 


(5) 
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v'(x)  -  ou(x)v(x),  v(0)  -  1. 

Observe  that  u(x)  satisfies  a  nonlinear  differential 
equation  whose  solution  is  determined  by  an  initial  condition 
as  compared  to  the  linear  equations  for  uR(y)  and  u^(y), 
determined  by  a  two-point  condition. 

5«  Connection  Between  the  Two  Approaches 

It  is  clear  that  by  suitable  choice  of  y,  we  can  obtain 
the  functions  u(x)  and  v(x)  from  the  functions  uR(y)  and 
uL(y).  Thus,  if  we  make  the  dependence  upon  x  explicit,  we 
have 

(1)  uR(y)  -  uR(y;x), 
uL(y)  -  uL(y;x), 

and 

(2)  u(x)  -  uR(x;x), 
v(x)  -  uL(0;x). 

Can  we,  however,  derive  the  internal  fluxes  uR(y)  and  vR(y), 
given  the  functions  u(x)  and  v(x)? 

To  accomplish  this,  we  combine  both  viewpoints.  Consider 
the  following  figure. 


0 


x 


Fig.  5 
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To  obtain  a  relation  between  uR(y),  u^(y),  and  u(y), 
we  consider  a  source  of  strength  u^(y)  per  unit  time  at  y. 
Then  the  steady-state  relation  is  clearly 

(?)  uR(y)  -  uL(y)u(y) . 

Similarly, 


(^)  u^(y)  -  v(x  -  y)  +  u(x  -  y)uR(y) . 


Hence , 

(5) 


solving  for  uR(y)  and  uL(y), 
u„(y)  - 

1  -  u(y)u( x  -  y) 

u  (y) - a - 

1  -  u(y)u(x  -  y) 


we  have 


It  follows  that  we  can  consider  u(x)  and  v(x)  as 
fundamental  functions  from  which  all  other  functions  can  be 
derived. 


6.  Semigroup  Properties 

Let  us  now  obtain  general  relations  connecting  u(x)  and 
v(x)  with  u(y),  v(y)  and  u(x  —  y),  v(x  —  y) .  The 
differential  equations  of  §4  are  particular  cases  of  these 
relations. 

Referring  to  the  figure  in  §5  and  tracing  the  multiply 
reflected  and  transmitted  fluxes,  we  sea  that 


For  a  definition  and  discussion  of  semigroups  see  1 


(1  )  u(x)  =  u(x  -  y)  +  v(x  -  y)u(y)v(x  -  y) 

+  v(x  -  y)u(y)u(x  -  y)u(y)v(x  -  y)  + 


p-iesB 
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Two  values  of  particular  interest  arc  y  *  A  and 
y  =  x  —  A.  The  value  y  -  x  —  A  leads,  as  A  — ►  0,  to  the 
foregoing  dJ  fferer.tial  equations,  and  the  value  y  =  A  to 

Stokes'  relations  5l],  a  matter  we  shall  discuss  again  beiow. 

These  results  show  that  we  can  replace  the  solution  of 
differential  equations  by  the  iteration  of  simple  transforma¬ 
tions,  Consequently,  these  relations  may  be  better  suited  for 
computational  purposes  than  the  foregoing  differential 
equations. 


7.  A  More  General  Imbedding 

The  foregoing  results  have  essentially  been  consequences 
of  the  observation  that  the  Internal  fluxes  uR(y)  and  u^(y) 
are  functions  not  only  of  y,  the  position  at  which  they  are 
measured,  but  also  of  x,  the  length  of  the  rod.  Hence,  we 
should  write f  as  already  noted  in 

(1)  uR(y)  -  u^yjx). 


uL(y)  -  uL(y;x) 
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Consider  now  the  more  general  situation  in  which  we  measure 
the  fluxes  at  y  due  to  a  source  at  an  internal  point  z. 

Source  Observer 

t  ? 

i - 1 - 1 - 1 

0  z  y  x 

Fig.  4 

The  right-hand  flux  at  y  should  now  be  denoted  by  uR(x,y,z) 
and  the  left-hand  flux  by  u^(x,y,z).  We  are  now  at  liberty 
to  allow  x,  y  and  z  to  vary,  either  independently  one  at 
a  time,  or  two  at  a  time,  or  all  three  together. 

We  see  then  that  there  are  at  least  three  different  ways 
in  which  we  can  imbed  a  particular  process  within  a  family  of 
processes.  Two  of  these,  variation  of  y  and  z,  lead  to 
linear  equations  with  two-point  boundary  conditions,  while  the 
third,  variation  with  respect  to  x,  leads  to  a  nonlinear 
equation  with  an  initial  value  condition.  Each  has  certain 
analytic  and  computational  advantages.  In  any  particular 
situation,  we  employ  the  formulation  which  is  most  convenient. 
For  further  detail,  see  [23] • 

8.  Energy  Dependence  (Wultigroup  Theory)  \j>j) 

Let  us  now  turn  our  attention  to  a  more  realistic  mathe¬ 
matical  model  in  which  we  assume  that  a  neutron  is  characteri¬ 
zed  by  an  energy  level  as  well  as  a  direction.  In  so  doing, 
we  have  our  choice  of  either  a  continuous  range  of  energies, 
or  a  finite  set  of  discrete  levels. 
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In  [23]  ,  we  have  di a cussed  the  continuous  version.  Let 
us  concentrate  upon  the  discrete  version  here,  since  this  is  a 
c Lee  of  greater  importance  from  the  computational  point  of 
view.  We  shall  begin,  as  before,  with  the  conventional 
formulation. 

The  internal  flux  at  y  is  now  described  by  two  vectors 


\lr) 

v1(y) 

(1) 

u(y)  - 

u2(y) 

• 

,  v(y)  - 

v2(y) 

• 

• 

• 

Vy). 

• 

• 

,Vy), 

Here  H  is  the  number  of  distinct  energy  levels  or  groups, 
u^y)  represents  the  flux  of  neutrons  in  the  i— th  level  to  the 
right,  and  v^(y)  the  corresponding  flux  to  the  left. 

Generalizing  the  foregoing  model  of  a  neutron  transport 
process,  we  suppose  that  various  interactions  such  as 
absorption,  fission  and  non-fission  collisions  and  so  on, 
result  in  neutrons  at  one  energy  level  being  transformed  into 
neutrons  at  other  levels. 

We  Introduce  four  matrices 

(2)  A  "  (ai j ) *  B"(bij)'  c-(cij)y  D  -  (dij) 

where 

(3)  aijA  •  the  expected  Incremental  number  of  neutrons 

at  the  i-th  level  in  the  right-hand  flux  at 
y  +  A,  per  neutron  at  the  J-th  level  in  the 
right-hand  flux  at  y. 
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to  within  terms  of  order  magnitude  o(a).  Similarly,  b.  .A 

-L  J 

denotes  the  incremental  contribution  from  left-hand  to  right- 
hand  flux,  from  right-hand  flux  to  left-hand  flux,  and 

d^jA  from  left-hand  flux  to  left-hand  flux.  It  should  be 
noted  that  in  the  completely  isotropic  case  the  matrices  A,  B, 
C,  D  are  closely  related. 

The  usual  conservation  considerations  lead  to  the  equations 


(4) 


u^y  +  A)  -  Uj_(y) 
V1  (y )  “  vt(y  +  A) 


N 


N 


=  A 


k  1  aiJUj(y)  +  A  1  biJVj(y)  +  °(A)» 

J  —  1  J  1 

N  N 

=  A  I  c,,u.(y)  +42  d.  .  v  (y )  +  o(a). 


J  =  1 


1J  J 


jrru  j 


for  1  ■  1,2, . . . ,H. 

Letting  h  — ♦  0,  we  obtain  the  following  vector-matrix 


equations  t 


(5)  !jH-Au  +  Bv, 

-  37  ”  011  +  Dv< 

for  0  £  y  £  x* 

As  before,  let  us  suppose  that  no  neutrons  are  incident 
at  0,  and  there  is  a  flux  of  intensity  b^  per  unit  time  of 
neutrons  in  the  i— th  level  at  x.  We  thus  obtain  the  two- 
point  boundary  conditions 


(6)  u(0)  -  0,  v(x)  *  b. 


If  the  rod  is  inhomogeneous ,  the  matrices  A,  B,  C  and 
J  will  depend  upon  y.  Although  there  is  no  difference  as 
far  as  the  functional  equation  technique  of  invariant  imbedding 
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la  concerned  between  the  treatment  of  the  homogeneous  and  In¬ 
homogeneous  ,  the  classical  treatment  Is  simplified  by  the 
assumption  of  constancy  of  A,  B,  C  and  D.  The  discussion 
below  applies  equally  to  constant  or  variable  matrices. 

Let  W(y)  be  the  matrix  solution  of 


(7) 


dV 

Hy 


V(0) 


I. 


To  solve  (5)  subject  to  (6),  we  suppose  that  v(0)  has  the  as 
yet  unknown  value  c.  Then,  the  solution  of  (5)  can  be 
written 


/  'v 


V  J 


Write 


(9) 


v(y)  - 


wu(y)  w12(y) 
W21(y)  W,j,(y) 


'22' 


where  each  is  an  N  x  N  matrix.  Using  the  terminal 

condition  v(x)  ■  b,  we  obtain  the  equation 


(10)  W22(x)o  -  b, 

which  determines  the  unknown  vector  c. 


9.  Computational  Aspects 

The  determination  of  c  in  (8.10)  requires  the  solution 
of  a  system  of  N  linear  equations  in  N  unknowns.  In 


p-ifc5b 


addition,  we  must  determine  the  N  x  N  matrix  W  using  the 
linear  differential  equation  in  (?)•  Fortunately,  since  the 
equation  is  linear,  we  can  determine  W(x)  one  column  at  a 

? 

time.  Hence,  instead  of  the  simultaneous  determination  of  N 
functions,  we  can  perform  N  determinations  of  N  functions. 

10.  Reflection  and  Transmission  Matrices 

Let  us  now  consider  the  foregoing  process  using  invari¬ 
ant  imbedding  techniques.  To  that  end  we  introduce  the  matrix 
R(x)  ■  (rj j(x) ) ,  where 

(1)  r.  .(x)  =  the  expected  flux  of  neutrons  in  state  i 

■*“  J 

reflected  per  unit  time  from  a  rod  of 
length  x  due  to  an  incident  flux  at  x  of 
unit  intensity  per  unit  time  in  state  J. 

The  same  type  of  reasoning  employed  in  the  one-dimensional  case 
yields  the  matrix  equation 

(2)  R(x  +  A )  *=  BA  +  ( I  +  AA)R(x)(l  +  DA)  +  R(x)CR(x)A  +  o 
In  the  limit  this  yields  the  Rlccati  matrix  equation 

(3)  R ' ( x )  =  B  +  AR  +  RD  +  RCR, 

with  the  initial  condition  R(0)  *  0. 

In  a  similar  fashion.  If  we  introduce  the  transmission 
matrix  T(x)  =  (t^j(x)),  where 
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(4)  t^(x)  =  the  exPected  flux  of  neutrons  In  state  1 

transmitted  per  unit  time  through  a  rod  of 
length  x  due  to  an  incident  flux  at  x 
of  unit  intensity  per  unit  time  in  state  J. 

Then,  we  obtain  as  before 

(5)  T ' ( x )  =  t(D  +  CR). 

11.  Computational  Aspects 

The  determination  of  R(x),  Dy  way  of  (10. 3),  requires 

2 

the  simultaneous  integration  of  N  nonlinear  equations  with 
the  initial  value  R(0)  =  0.  This  is  a  far  more  complicated 
operation  than  that  of  solving  N  sets  of  N  linear  equations, 
but,  in  recompense,  it  avoids  the  task  of  solving  N  simul¬ 
taneous  linear  equations. 

Furthermore,  let  us  note  that  once  R(x)  has  been 
determined,  we  have  resolved  the  transport  process,  determina¬ 
tion  of  internal  and  external  fluxes,  for  a  set  of  rods  of 
increasing  length.  On  the  other  hand,  the  conventional 
method  based  upon  linear  equations  yields  the  solution  for  one 
length  at  a  time. 

12.  Criticality 

Let  us  turn  to  a  discussion  of  one  of  the  most  important 
phenomena  associated  with  neutron  transport  and  multiplication, 
namely  criticality.  As  the  length  of  the  rod  Increases,  the 
intensity  of  internal  and  emergent  flux  increases  and  becomes 
infinite  as  a  certain  critical  length  is  attained. 
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To  determine  the  critical  length  for  the  energy-independent 
case,  let  us  begin  with  the  linear  equations  of  §2.  Eliminating 
uL(y),  we  obtain  the  equation 

U)  UR  (y)  =  =  “  a2uR(y)* 

We  take  a  constant  for  simplicity.  The  general  solution  Is 

(2)  uR(y)  =  ci  8ln  °y  +  c2  G0S  oy* 

Using  the  two-point  boundary  conditions  of  (2.5  ) $  we  readily 
obtain  the  equation 

(3)  uR( y )  *  sin  oy/cos  ox. 

We  see  then  that  uR(y)  and  uL(y)  are  infinite  for 
0  <  y  <  x  when  x  =  w/ 2a.  This  is  the  critical  length  for  the 
simple  neutron  multiplication  process  we  have  set  up. 

Turning  to  the  equation  for*  the  reflected  flux  obtained 
via  Invariant  Imbedding,  we  have 

(4)  u'(x)  =  o(l  +  u2),  u(0)  =  0, 

whence 

(5)  u(x)  =  tan  ox. 

Once  again,  we  see  that  x  «  ir/2a  is  the  critical  length. 

As  pointed  out  by  D.  McOarvey,  we  can  use  (6.1)  or  (6.2) 
to  obtain  the  critical  length.  In  place  of  asking  for  the 
value  of  x  which  makes  u(x)  Infinite,  it  is  sufficient  to 
ask  for  the  value  of  x  which  makes  u(x)  »  1,  and  then 
double  x. 
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lj.  Critical! ty--Multlgroup  Case 

It  is  in  the  determination  of  critical  length  in  the 
energy-dependent  case  that  the  classical  formulation  encounters 
real  trouble.  To  find  the  value  of  x  which  yields  infinite 
flux,,  we  must  solve  the  determinantal  equation 

(1)  det  (W22(x))  -  0. 

Let  us  note,  once  and  for  all,  that  when  we  speak  of  the 
critical  value,  we  mean  the  smallest  value  of  x  which  yields  an 
infinite  flux.  Prom  the  physical  point  of  view,  the  problem 
of  determining  expected  fluxes  is  meaningless  when  the  length 
of  the  rod  exceeds  the  critical  value.  The  higher  values  of 
x  which  yieid  Infinite  values  are  connected  with  the  higher 
characteristic  values  associated  with  the  two-point  boundary- 
value  problem.  They  do  not  appear  to  have  any  physical 
significance,  although  this  Is  always  a  dangerous  statement. 

On  the  other  hand,  when  we  go  over  to  a  more  sophisticated 
discussion  concerning  probabilities  of  fluxes  of  various 
intensity,  and  probability  of  fission,  then  it  becomes  quite 
significant  to  consider  rods  of  greater  than  critical  length. 

There  are  a  number  of  interesting  mathematical  problems  In  this 
area  which  have  oeen  considered  in  detail  by  McGarvey,  m. 
and  Mullikln  and  Snow,  [25]  .  We  shall  disouas  than  briefly  below. 

Returning  to  the  equation  In  (l),  we  see  that  If  N,  the 
number  of  groups,  is  of  any  size,  say  10  or  20,  the 
problem  Is  not  trifling.  If  N  »  50  or  100,  we  cannot  consider 
a  solution  along  the  foregoing  lines  to  be  satisfactory,  for 
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a  number  of  reasons  which  are  familiar  to  numerical  analysts. 

The  Invariant  imbedding  technique  requires  the  integration 

2 

of  N  simultaneous  differential  equations  which  are  quadratl- 
cally  nonlinear.  This  integration  is  pursued  until  some 
element  in  the  matrix  R(x)  becomes  infinite.  To  begin  with, 
let  us  discuss  the  dimensional  aspects.  A  computation  of  this 
type  for  N  «=  10  or  20  is  completely  routine  for  modem 
digital  computers,  and  one  of  this  nature  for  N  «  50  is 
large,  but  feasible.  For  the  machines  that  will  be  operational 
within  a  few  years,  values  of  N  such  as  100  or  200  will 
be  routine. 

Now  let  us  turn  to  the  integration  of  the  differential 
equations  until  a  singularity  occurs.  Clearly  this  is  not  a 
routine  operation  if  accuracy  is  desired.  There  are  several 
things  that  we  can  do.  First  of  all,  we  can  observe  that  as 
x  approaches  xQ,  the  critical  value,  we  have  an  asymptotic 
behavior  of  the  form 


(2) 


rU(x) 


T- -ii 

(x  -  X 


0 


where 

(3) 


'u 


1 0. 

1 


and  some  ^  >  0. 

(x  -  xQ) 

~ - 

8u 


Hence 


when  s^  >  0.  This  linear  behavior  can  be  used  to  predict 
the  value  of  with  great  accuracy.  Furthermore,  the  fact 

that  there  are  functions  i^j(x)  will  enable  us  to  deter¬ 

mine  Xq  with  even  greater  accuracy. 


P-1858 

29 


Secondly,  we  can  use  McGarvey's  observation,  pointed  out 
In  the  section  on  criticality  for  the  eimple  energy- independent 
case.  In  place  of  finding  the  first  value  of  x  for  which 
R(x)  is  singular,  we  can  ask  for  the  first  value  of  x  for 
which  the  matrix  R(x)  has  its  largest  characteristic  root 
equal  to  one.  If  this  value  is  ,  the  critical  value  will 
be  2x^ . 

Since  the  matrix  H(x)  is  a  positive  matrix,  or  at  least, 
nonnegative,  we  know  that  there  will  be  one  root  of  largest 
absolute  value  which  is  real  [5^]*  By  slight  perturbation  of  the 
transition  matrices  A,  B,  C,  and  D  we  can  actually  ensure 
that  all  the  entries  in  R(x)  are  positive,  which  means  that 
the  root  with  largest  absolute  value  will  actually  be  positive. 

There  are  now  available  a  number  of  simple  and  efficient 
techniques  for  determining  thl3  root,  the  Perron  root,  of  a 
positive  matrix.  Furthermore,  since  clearly  R(x)  has  mono- 
tonically  Increasing  elements,  we  can  use  various  interpolation 
methods  to  locate  the  position  of  this  root  very  accurately. 

A  large  number  of  questions  in  the  theory  of  branching  processes 
can  be  reduced  to  the  problem  of  determining  the  largest 
characteristic  roots  of  positive  operators;  see  Bellman-Harris, 

[26]  ,  Blrkhoff ,  [27]  . 

In  any  case,  this  method  seems  far  superior  to  that  of 
finding  the  roots  of  a  determinants 1  equation  of  high  degree. 

It  would  seem  that  invariant  imbedding  techniques  have  a 
distinct  advantage  as  far  as  the  determination  of  critical 
parameters  is  concerned. 
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14.  Extrapolation  over  Multigroups 

One  way  of  determining  the  critical  length  with  great 
accuracy  la  based  upon  the  use  of  a  large  number  of  energy 
levels.  It  1 8  reasonable  to  suspect  that  closer  and  closer 
values  to  the  true  value  will  be  derived  as  we  use  finer  and 
finer  subdivisions  of  the  energy  range.  Consequently,  we  can 
use  the  following  extrapolation  method.  Solve  the  problem  for 
N  »  10,  for  N  =  20,  N  ■  30,  and  so  on,  until  we  reach  the 
limits  of  the  computer.  Using  the  successive  values  obtained 
for  the  critical  length,  we  can  extrapolate  to  N  =  ®,  and 
thereby  obtain  a  more  precise  value. 

Here  is  where  an  analysis  of  the  precise  asymptotic  form 
as  N  — ►  00  will  be  very  valuable.  With  the  aid  of  an 
analytic  representation  of  the  critical  length  as  a  funotion 
of  N,  we  ean  use  superior  extrapolation  procedures. 

Of  course,  it  is  seldom  possible  to  obtain  cross  sections 
and  other  physical  parameters,  as  continuous  functions.  Hence 
the  limiting  case,  N  — *  00  ,  is  often  of  greater  mathematical 
Interest  than  it  is  of  physical  Importance. 

15.  Multidimensional  Transport  Theory- -Slab  Case 

Leaving  the  physically  cramped  but  mathematically  com¬ 
fortable  confines  of  the  one-dimensional,  let  us  begin  our 
investigation  of  the  more  significant  multidimensional  pro¬ 
cesses  oy  considering  a  neutron  transport  process  taking  place 
in  an  infinite  slab  contained  between  the  planes  y  =  0  and 

y  =  x  in  three  space.  As  usual,  surrounding  the  slab  is  a 
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vacuum  which  inean3  that  a  neutron  leaving  the  slab  at  either 
boundary  never  returns. 


Pig.  5 


A  classical  formulation  of  this  problem  leads  in  the  iso¬ 
tropic  case  to  the  equation 


(1) 


i  2£  + 

c  at  r  a y 


of 


■  ¥/; 


f ( JCyx’  ,t)d/4' 


where  c  is  the  constant  neutron  velocity,  o  is  the  constant 
collision  cross-section  and  k  is  the  average  number  of 
neutrons  emerging  from  a  collision.  As  usual,  u  is  the 
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cosine  of  the  angle  between  the  direction  of  motion  of  the 
particle  and  the  positive  y-direction,  and  f(y,/^t)  is  the 
density  of  neutrons  at  y  traveling  in  direction  ^u.  at 
time  t. 

There  are  boundary  conditions  at  y  =  0  and  y  *  x, 
arising  from  the  fact  that  particles  may  not  reenter  the  slab 
once  they  have  emerged.  In  the  steady-state  situation,  of  the 
type  we  have  so  far  been  considering,  (l)  takes  the  form 

(2)  />!£  +  or  -  y*1  . 

A  rigorous  treatment  of  Equations  (l)  and  (2)  requires  deep 
analysis;  see  [20,2lj  . 

Let  us  consider  this  problem  using  invariant  imbedding 
techniques.  To  simplify  our  initial  presentation,  let  us 
return  to  neutrons  which  are  independent  of  energy,  but  do, 
however,  possess  directions  of  motion.  Assume,  as  indicated 
in  the  picture  above  that  there  is  a  plane  parallel  flux  in 
direction  0  per  unit  area  per  unit  time  incident  at  x, 
and  that  we  are  given  the  various  probabilities  of  absorption, 
scattering  and  fission  collisions,  and  the  resultant  angular 
distribution  of  neutrons. 

The  type  of  reasoning  used  in  the  previous  sections 
enables  us  to  derive  a  functional  equation  of  the  form 

(5)  &-*»>. 

where  T  is  a  quadratic  operation,  for  the  function 
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(4)  u(x,e,^)  -  the  reflected  flux  per  unit  area  on  the 

surface  in  the  ^-direction  per  unit  time 
as  a  result  of  a  unit  incident  flux  per 

• 

unit  area  on  the  surface  per  unit  time. 

There  la  no  need  for  us  to  go  into  the  details  for  three 
reasons.  In  the  first  place,  we  shall  derive  similar  equations 
below  for  cylindrical  and  spherical  geometries.  The  second 
reason  we  shall  discuss  immediately  below.  Finally,  we  shall 
discuss  this  problem  from  another  physical  viewpoint  in  Part  V. 

16.  Equivalence  of  One-dimensional  Energy- dependent  Case  and 

Angular-dependent,  Energy- Independent  Slab  Case 

What  is  Important  is  the  observation  that  the  transport 
process  for  a  one-dimensional  rod  with  energy-dependence  is 
abstractly  equivalent  to  the  process  for  the  Blab  with  discrete 
angular  dependence,  but  no  energy  dependence. 

In  both  cases,  we  have  a  finite  number  of  "states”  and 
mechanisms  for  transforming  a  neutron  from  one  state  to 
another.  Another  advantage  of  this  formulation  lies  in  the  fact 
that  the  inclusion  of  energy  dependence  in  the  slab  merely 
increases  the  number  of  states,  without  at  all  changing  the 
mathematical  formulation. 

17.  Cylindrical  Regions 

The  infinite  slab  is  stratified  by  considering  it  to  be 
composed  of  a  series  of  strata  of  which  the  stratum  between 
y  =  x  and  y  =  x  —  A  is  typical  (Figure  6). 

f~~  — • 

Throughout  this  paper  we  measure  fluxes  with  respect  to  the 
geometrical  areas  on  which  they  impinge,  rather  than  with  respect 
to  a  plane  normal  to  the  beam.  For  a  discussion,  see  §59* 
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Fig.  6 

In  analogous  fashion,  we  can  stratify  other  regions  with 
various  types  of  symmetries.  Consider,  for  our  first  example 
of  this,  an  infinite  cylindrical  region,  whose  cross-section 
d(r)  is  dependent  only  upon  the  radial  coordinate  r.  Let 
us  suppose  that  neutron  production  is  energy  independent  and 
isotropic,  with  k  neutrons  emerging  after  each  collision. 

Given  an  Incident  flux  of  one  neutron  per  unit  area  per 
unit  time  on  the  surface  at  angles  (6,^),  we  wish  to  deter¬ 
mine  the  reflected  flux  ) .  As  usual, 

yu.  m  cos  0  (Figure  7). 

The  imbedding  is  now  performed  by  considering  the 
cylindrical  region  to  be  composed  of  a  sequence  of  infinitesi¬ 
mal  cylindrical  shells.  In  cross-section,  they  appear  as  in 
Figure  8. 
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(For  ingoing  particles ,  0  is  measured  with  respect  to  the 

inward  normal.) 


Fig.  7 


Fig.  8 


Referring  to  the  foregoing  two  figures,  and  adding  up  effects 
as  before,  we  obtain  the  functional  equation  (see  [47]) 
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no  ,t')  =  0. 


To  compute  d^/dr  we  must  note  that  /^  Itself  Is  really 
a  function  of  r,  and  the  same  i3  true  of  /l '  .  Upon  taking 
this  Into  consideration  we  find 

to\  d^  dv  dv  1  -/^2  1  -  A'2 

dr  dr  /or  d/o'  /o'r 

We  shall  discuss  the  corresponding  result  for  spherical 
regions,  and  then  discuss  the  computational  significance  of 
these  results. 

l8.  Spherical  Regions 

As  our  nexc  example,  consider  a  sphere  composed  of  trans¬ 
port  material  whose  cross-section  is  dependent  upon  the  radial 
coordinate  p  alone.  As  indicated  in  the  following  figure,  we 
introduce  an  angular  coordinate  a,  cos  a  ■  v,  and  suppose 
that  we  have  a  conical  flux  of  neutrons,  with  direction  v, 
incident  unifonnly  over  the  surface  of  the  sphere,  one 
neutron  per  unit  area  per  unit  time.  We  wish  to  determine  the 
reflected  flux  in  direction  v’,  ^f(p,v,v'). 


P-1858 

57 


if 

:  / 

i  / 

p  A 

/  V 

Pig.  9 

The  usual  analysis  yields  the  equation 

(1)  ip  +  M?p~~  *  ~v  V  ^ 

-  +  a’^£i-c/71  V(p»v",v')dv" 

-  d(p)(i  +  7r)f(p.*,V) 

+  |  dfp)^1  dv"  (l  +  f(p,v’",v')dv'"j. 

19.  Critical  Hass 

The  critical  reader  nay  seriously  question  the  value  of 
the  results  obtained  in  the  two  previous  sections,  since 
Incident  fluxes  of  the  type  we  have  employed  are  seldom  found. 

This  Is  certainly  a  valid  criticism. 

There  Is,  however,  one  quite  Important  case  In  which  we 
can  profitably  use  this  type  of  flux,  and.  Indeed,  whatever 
type  of  flux  Is  most  convenient.  This  is  the  determination  of 
critical  mass.  Zt  Is  possible  to  convince  oneself  that  whatever 
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radius  Is  critical  for  one  type  of  flux  will  be  critical  for 
any  other  type  of  steady-state  flux. 

20,  More  Oeneral  Fluxes 

The  aajne  persevering  reader  may  also  ask  why  we  have  not 
used  Invariance  principles  directly  for  more  general  flues. 

This  can  be  done.  What  has  held  ua  back  has  been  dimensionality 
difficulties.  Consider,  for  example,  a  two-dimensional  slab 
In  which  we  consider  an  Incident  flux  of  un< t  intensity  per 
unit  time  at  an  angle  9  at  a  particular  point,  say  z  =  0. 


Fig.  10 

We  then  introduce  the  flux,  u(6,P,z,x),  a3  the  reflected 
flux  at  angle  p*  per  unit  time  at  the  point  a  distance  z 
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from  the  point  of  Incidence.  We  then  obtain  the  same  type  of 
equation  for  u  as  before,  with  the  difference  that  u  now 
depends  upon  one  additional  variable  z.  This  Increase  in 
dimensionality  introduces  formidable  computational  difficulties 
due  to  the  enormously  increased  memory  requirements. 

From  the  analytic  point  of  view,  there  is  no  difficulty 
in  considering  realistic  fluxes.  From  the  computational  point 
of  view,  these  more  realistic  problems  require  new  techniques 
and  bigger  and  faster  machines.  For  a  possible  line  of 
approach,  see  [28]  . 

21.  Vol terra  verous  Fredholm  Equations 

The  equations  we  have  obtained  in  the  foregoing  sections 
using  invariant  imbedding  techniques  have  invariably  been 
nonlinear,  as  compared  to  the  linear  equations  of  classical 
transport  theory.  Considering  the  fact  that  linear  equations 
with  their  superlative  superposition  properties  are  difficult 
enough  to  analyze  theoretically  and  resolve  numerically,  why 
do  we  wish  to  introduce  nonlinear  equations?  Although  we  have 
gnawed  around  the  edges  of  this  question  in  previous  pages,  let 
us  now  make  it  the  principal  course. 

Our  answer  can  be  put  in  very  simple  terms  1  We  wish 
uniformly  to  replace  two— point  boundary  value  problems,  and 
boundary-value  problems  in  general,  by  initial  value  problems. 
This  is  equivalent  to  replacing  Fredholm— type  integral 
equations  by  Yolterra— type  integral  equations.  Naturally,  the 
equations  will  be  in  different  variables. 
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This  is  not  a  new  idea,  and,  indeed,  la  one  that  has  been 
proposed  before,  and  used  with  some  success.  The  approach  we 
use,  however,  based  upon  invariance  principles,  and  extending 
that  of  Aabarzumian ,  [29]  ,  and  Chandrasekhar,  [}0]  ,  is  ^uite 
different  from  these  mentioned,  and  quite  unlike  any  earlier 
methods . 

We  are  engaged  in  this  program  for  two  reasons.  In  the 
first  place,  Volterra— type  equations  lead  to  Iterative 
algorithms  which  are  simpler  for  digital  computers  than 
algorithms  based  upon  the  solution  of  linear  systems  of 
equations.  Secondly,  and  the  two  are  intimately  related, 
functional  equations  of  the  type  we  derive,  despite  their 
nonlinearity,  are  easier  to  treat  as  far  as  existence  and 
uniqueness  are  concerned.  These  la3t  topics,  however,  we 
have  bypassed  here  since  we  are  at  the  moment  principally 
interested  in  exhibiting  the  mechanics  of  the  formulation  of 
the  classical  particle  processes  in  these  new  terms. 

22,  Stokes  Relations 

It  is  interesting  to  find  that  the  reflection  and  trans¬ 
mission  matrices,  and,  more  generally,  the  reflection  and 
transmission  functions,  are  related  to  each  other  in  algebraic 
fashion.  We  already  have  the  Riccati  differential  equation 
for  R(x)  and  the  differential  equation  for  T(x)  in  which 
R(x)  enters  (see  §10) . 


We  shall  refer  to  these  new  relations  as  Stokes ' 
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relations  since  the  first  Identities  of  this  type  connecting 
reflection  and  transmission  coefficients  were  discovered  by 
Stokes  in  work  on  light  rays  impinging  on  slab*. 

In  §  6  ,  we  discussed  the  semigroup  properties  of  the 
functions  u(x)  and  v(x)  obtained  for  the  3imple  one- 
dimensional  energy— independent  process.  Recognition  of  these 
transformation  properties  is  also  found  in  Redheffer,  QjlJ  . 
and  for  the  matrix  case  as  well.  To  obtain  the  following 
relations,  we  use  the  special  case  of  these  relations  in  which 
the  stratification  instead  of  being  as  before, 

i - 1 - 1 

0  x  x+A 

Pig.  11 


is,  instead, 

I - f 

0  A 

Pig.  12 

where  A  is  an  infinitesimal. 

The  usual  counting  process  yields  the  relation 

(1)  rt(x  +  A)  =  R(x)  +  T(x)BT(x)A  +  o(a). 


H 

X-'  A 


or 
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(2)  R ' (x)  -  TBT. 

In  view  of  Equation  (10.3)*  we  conclude  that 
(5)  TBT  "  B  +  AR  +  RD  +  RCR, 

a  relation  which  is  not  too  easy  to  verify  directly. 
Similarly,  we  find  that 

(4)  T(x  +  A )  *=  ( I  +  DA)T(x)  +  R(x)BT(x)  A  +  o(A), 
whence 

(5)  T'(x)  -  (D  +  RB)T. 

This  in  conjunction  with  Equation  (10.5)  yields  the  further 
result 

(6)  T(D  +  CR)  -  (D  +  RB)T. 

23.  Probabilities 

Por  various  purposes,  it  is  desirable  to  have  a  more 
precise  picture  of  the  transmitted  and  reflected  flux  than 
that  furnished  by  the  expected  values.  This  ia  particularly 
the  case  when  x  is  Just  less  than  critical,  when  it  assumes 
the  critical  value,  and  when  x  is  greater  than  critical, 
the  case  of  supercritical! ty. 

Let  us  then  talk  about  the  probabilities  of  events, 
rather  than  the  expected  events.  In  place  of  a  steady-etate 
situation,  let  us  suppose  that  one  neutron  is  incident  at  x. 
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the  right  end  of  our  one-dimensional  rod,  at  time  zero.  We 
shall  call  this  neutron  the  trigger  neutron. 

Let  us  define  the  set  of  probabilities 

(1)  p^(x)  "  Prot5abillty  that  n  neutrons  are 

reflected  from  a  rod  of  length  x  over 
all  time  as  a  result  of  one  trigger 
neutron  incident  at  x  at  time  zero. 

In  order  to  obtain  a  set  of  differential  equations  for  these 
functions,  we  observe  that  the  trigger  neutron  incident  upon 
a  rod  of  length  x  +  A  either  has  a  fission  collision  in  the 
initial  length  [x  +  A,x]  or  it  does  not.  If  it  does  not, 
consider  the  neutrons  that  emerge  from  the  length  [0,xj.  If 
k  of  these  emerge  over  all  time,  then  at  most  one  of  these, 
to  terms  in  o(a),  can  have  a  fission  collision  in  the 
Interval  [x,x  +  A^  .  Consideration  of  these  possibilities, 
together  with  that  of  an  initial  fission  collision,  leads  to 
the  relation 

(2)  Pn(x  ♦  A)  «  (1  -  OA)jpn(x)(l  -  nOA)  + 

♦  +  °(  a  ) , 

n  ■  0,1,2,...;  P_^('*)  "  0. 

Passage  to  the  limit  yields  the  infinite  system  of 
differential  equations 
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(3)  Pn 1 (x)  ■-(**♦  l)^Pn(x)  +  dp^  +  dZ  kpk(x)pn_k(x)# 

k*l 


n  ■  0|1|2, • • • )  p_^ ( x )  ■  0 • 

Physical  considerations  yield  the  obvious  boundary  conditions 


pn(°) 


1,  n  -  0 
0,  n  yi  0. 


In  order  to  test  the  feasibility  of  the  computational 
solution  ‘of  systems  of  this  type,  and  to  obtain  some  idea  of 
the  time  required,  the  first  forty  of  these  equations  were 
solved  numerically  using  a  rather  old-fashioned  machine,  the 
RAND  Johnniac.  Graphs  of  the  results  may  be  found  in  [ ] 
and  [33]. 

A  number  of  interesting  questions  concerning  the  analytic 
behavior  of  the  Pn(x)  as  x  approaches  the  critical  value 
arise  as  a  result  of  these  computations. 

If  we  introduce  the  generating  function 

00 

(5)  u(x,s)  -  Z  p  (x)sn, 

n-0  n 


we  readily  obtain  the  quasllinear  partial  differential 
equation 

(6)  ux  -  d(s  -  l)u  +  ds (u  -  l)u3, 

with  the  initial  condition 


(7) 


u(0,s)  ■  1. 
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This  equation  has  been  studied  in  detail  by  Hullikin  and  Snow, 
[25j  •  The  solution  exhibits  a  very  strange  behavior  as 
x  goes  through  the  critical  value.  As  we  shall  see  below, 
the  equation  for  the  generating  function  can  be  obtained 
immediately  by  means  of  functional  equation  techniques. 

It  is  important  to  insert  a  word  of  caution  about  the  use 
of  these  techniques.  Methods  based  upon  expected  values  can 
be  used  with  a  carefree  abandon  as  long  as  the  length  of  the 
rod  is  less  than  critical.  On  the  other  hand,  the  probabilis¬ 
tic  equations  obtained  above  hold  for  any  value  of  the  rod. 
They  must  be  considered  the  basic  equations. 

As  long  as  x  is  less  than  the  critical  length,  we  have 
the  condition 

(8)  2  m  1 

n-0 

in  addition  to  the  initial  values  of  (4).  As  soon  as  x  is 
greater  than  the  critical  length,  this  equation  no  longer 
holds.  This  is  due  to  the  fact  that  there  is  now  a  new 
probability  to  be  added,  the  probability  of  an  infinite  flux. 
Lack  of  recognition  of  this  fact  can  lead  to  paradoxical 
conclusions.  A  thorough  discussion  of  this  phenomenon 
will  be  found  in  the  paper  by  Mullikin  and  Snow  cited 
above . 

Finally,  let  us  note  that  similar  equations  can  be  ob¬ 
tained. for  the  transmission  probabilities,  see  [32]  . 
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24,  Analogy  between  Critical  Length  and  Initiation  of  Shock  Wave 
The  equation  for  the  generating  function,  given  above  in 
^23,  la  analogoua  to  the  equation 

(1)  ut  -  uux, 

where  t  is  now  a  time  variable  and  ace  variable, 

used  by  Courant-Hilbert,  and  others,  as  sample  of  how  a 
discontinuity  can  arise  in  the  behavior  of  the  solution  of  a 
differential  equation.  This  type  of  discontinuity  is  similar 
to  that  which  is  observed  in  the  behavior  of  blast  waves  and 
is  called  a  "shock." 

It  is  interesting  then  to  observe  this  analogy  between 
the  onset  of  a  shock  as  a  function  of  time  and  the  onset  of 
criticality  as  the  length,  or  radius,  is  increased.  Analogies 
of  this  type  are  useful  since  knowledge  gained  in  one  area  can 
then  be  easily  transplanted  to  another. 

We  3hall  observe  a  further  analogy  subsequently.  Just  as 
the  presence  of  the  smallest  degree  of  viscosity  destroys  the 
pure  shock,  so  the  presence  of  the  slightest  neutron— neutron 
Interaction  can  destroy  criticality  (see  [3^]). 

23.  Description  of  a  generalized  Transport  Process 

As  we  shall  see  in  a  moment,  the  imbedding  technique, 
which  we  have  used  so  far  only  for  one-dimensional  rods, 
slabs,  cylindrical  and  spherical  regions,  can  be 
considerably  extended.  Let  us  now  formulate  a  trans¬ 
port  process  in  general  terms.  Let  a  family  of  surfaces  in 


n-diraensional  space  be  characterized  by  a  single  continuous 
parameter  7)  ^  0.  The  surface  corresponding  to  V  -  0  (it 
may  be  a  degenerate  surface)  will  be  considered  a  bounding 
surface,  and  the  parametrization  will  be  such  that  the  region 
included  between  7?-  0  and  7?*  ^  will  be  included  by  the 
region  between  0  and  ??  -  ^  >  7?  It  is  clearly  not 

necessary  to  consider  the  topological  properties  of  these 
surfaces  and  regions  in  detail.  Such  ideas  as  outward  and 
inward  directions,  area  of  a  surface,  and  so  on,  can  be  taken 
intuitively  as  far  as  our  purposes  are  concerned.  Ifae  family 
of  surfaces  will  be  assumed  to  partition  continuously  all  or 
part  of  the  n-dimensional  space  into  a  set  of  strata. 

Por  example,  the  family  of  surfaces  may  be  the  set  of  all 
spheres  centered  at  the  origin  in  three-dimensional  space. 

Here  7]  is  r,  the  radius  of  the  sphere,  and  0  is  a 

degenerate  bounding  surface,  a  point  sphere.  Again,  the 
family  may  be  the  set  of  all  vertical  lines,  in  the  two-dimen¬ 
sional  space,  to  the  right  of  the  vertical  axis.  In  this  case 
we  can  choose  7?  to  be  x,  and  x  *  0  Is  a  nondegenerate 
bounding  surface. 

By  a  "particle"  we  shall  understand  a  state  vector  S 
depending  on  the  parameter  7?.  The  state  vector  contains 
information  regarding  the  direction  of  motion,  energy,  specific 
location  on  77,  the  type  or  types  of  physical  particles  that 
we  are  discussing  (in  case  there  are  particles  other  than 
neutrons),  and  any  other  properties  which  we  may  choose  to 
Include  m  the  process. 
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The  stratum  (??,  7)+  A) ,  A  >  0,  will  be  assumed  to  con¬ 
tain  a  medium  which  permits  a  transport  process.  A  particle 
passing  through  the  stratum  may  engage  in  interactions  of 
both  a  deterministic  and  a  stochastic  nature.  As  mentioned 
above,  the  distinction  between  these  is  really  a  matter  of 
mathematical  convenience.  The  deterministic  interactions 
will  produce  effects  which  will  be  proportional  in  magnitude 
to  A  plus  a  term  o(A);  the  stochastic  interactions  will 
have  probability  of  occurrence  proportional  to  A  plus  o(a) 
within  the  stratum  (7?,/?+  A).  These  interactions  will  have 
no  effect  on  the  transport  medium,  but  will  result  in  a  trans¬ 
formation  on  the  state  vector  S,  in  some  cases  transforming 
it  into  two  or  more  such  vectors  (fission),  in  some  cases 
annihilating  It  (absorption),  in  other  cases  leaving  one 
vector  as  before.  Hence,  the  movement  of  a  particle  within 
the  medium  between  0  and  71  can  be  thought  of  as  a  sequence 
of  transformations  on  state  vectors,  together  with  the  creation 
and  annihilation  of  vectors.  Subsequently,  in  §36,  we  shall 
consider  a  process  in  which  the  medium  is  changing  as  the 
process  continues. 

Ve  shall  now  investigate  the  following  problem.  Let  a 
flux  of  "particles”  specified  by  S  impinge  on  ??.  What  are 
the  number  and  nature  of  the  particles  (state  vectors) 
emergent  from  7/f  Often  the  source  will  be  given  in 
particles  per  unit  time,  in  which  case  we  shall  seek  the 
number  of  particles  emerging  from  7J  per  unit  time.  In  some 
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instances,  we  shall  be  concerned  with  the  number  and  nature  of 
particles  emergent  from  /y  »  0  per  unit  time?  The  former 
particles  we  shall  call  "reflected,"  the  latter  we  call  "trans¬ 
mitted."  It  is  clear  that  in  some  cases,  such  as  for  the  sphere 
mentioned  above,  the  second  problem  is  ill-posed.  We  shall 
investigate  only  the  problem  of  reflection,  since  the  formulas 
for  transmission  can  be  obtained  in  a  similar  fashion. 

26.  The  Expected  Value  Equation 

We  begin  by  considering  not  the  state  vector  Itself,  but 
its  expected  value.  Consider  a  stream  of  particles,  one  per 
unit  area  per  unit  time,  in  state  S,  impinging  upon  the 
surface  7 ).  We  ask  for  the  expected  number,  or  flux,  of 
particles  in  state  S'  reflected  per  unit  area  per  unit  time. 
Denote  this  flux  by  Y(??, S,S').  Por  the  present  discussion  we 
shall  assume  that  S  and  S'  contain  no  information  about 
the  specific  location  on  7?.  It  is  evident  that  this 
assumption  imposes  a  strong  symmetry  requirement  upon  both 
the  surface  7?  and  the  Impinging  flux. 

Let  the  probability  that  a  particle  in  state  S  passing 
through  (7?*7?+  A)  suffers  a  collision— the  stochastic 
process — be  given  by  P(??,S)A  +  o(A).  Let  T(77,S,S’)  be  the 
average  number  of  particles  in  state  S'  resulting  from  this 
interaction  and  transmitted  to  y  +  A.  Let  R(^7,S,S')  be  the 
average  number  in  state  S'  reflected  back  to  77.  Then, 
proceeding  as  above,  we  find 
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(1) 


Y(t*-a,s,s') 


+ 

+ 

+ 


P(?,S)AR(7?,S,S') 

P(J!,S)4/  T(^S,S")t'(^.S",S’)<lS'' 

[l  -  P(7?,3)A]f(?7,3,S')[l  -  P(54S')A] 
[l  -  P(?,S)A]  dS'’Y(^3,S’')AP(^S") 

[*(y,SM,S')  +y  R(?.S\S''')Y(?,S"’ 


,S')dS"  ' 


♦  o(a)  . 


The  first  terra  on  the  right-hand  side  of  this  equation 
represents  the  contribution  to  the  reflected  flux  from  the 
particles  which  are  reflected  immediately  from  the  stratum  of 
thickness  A,  and  the  second  arises  from  those  which  pass 
through  this  stratum  but  change  state,  giving  reflected 
particles  from  [o  ,Pp]  which  then  pass  through  the  stratum 
without  Interaction.  The  third  term  is  produced  by  particles 
which  pass  unaffected  through  [/?,/?  +  A]  and  give  rise  to 
reflected  particles  from  [ 0,7J\  which  again  pass  through  the 
stratum  without  interaction.  The  last  term  accounts  for 
particles  which  enter  [0,?3  without  interaction  in  [/?,  V*  A], 
but  whose  reflected  flux  does  have  an  interaction  in  [??,??  +  A], 
Some  of  this  flux  is  transmitted  through  this  stratum,  the  rest 
is  returned  to  [0>2i  40(1  re-reflected. 

Here  all  the  states  on  the  right-hand  side  are  taken  at 
7J,  and  the  integration  over  states  is  symbolic .  The  trans¬ 
port  equation  resulting  when  A  — ►  0  is  then 
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<2>  P(7?,S)B(^,S,S')  +  dS"T(j),S,S")f(^,S",S-) 

s 

-  [P(7?.S)  +  P(^,S')]T(?.S,S’) 

dS"f(7?,3,3")P(^,S")|T(^S",S') 

+y?  dS1"R(7,S\S"'m9,S'",S 

It  may  be  shown  that  (2)  Includes  Equations  (17.1)  and 
(18.1)  as  special  cases.  By  choosing  different  geometries  and 
different  meanings  for  S  it  is  possible  to  write  down  a  great 
variety  of  particular  transport  equations  using  this  general 
result.  For  more  details,  see  [4]. 

27.  A  Simple  Stochastic  Case 

Before  attempting  to  write  down  a  general  stochastic 
functional  equation,  we  consider  a  simple  example,  the  first  of 
^2.  We  assume,  for  further  simplicity,  that  in  a  collision 
exactly  two  neutrons  emerge,  one  going  to  the  right,  one  to  the 
left.  We  let  one  neutron  enter  the  bar  at  time  zero. 

Let  ^U^(x)j,  1  »  1,2,3,...,  be  a  sequence  of  random 

variables  denoting  the  number  of  neutrons  reflected  from  the 
bar  in  all  time.  Let  j^F^(A)j,  i  -  1,2,3,...,  be  another 
sequence  of  random  variables  with 

/.x  (l  if  a  collision  occurs  in  a  segment  of  length  A 

F(1)(A)  - 

10  otherwise. 


0 


- 1- 

x—iA 

Fig.  15 


~! 

x  ♦—  (Initiating 
particle) 


P— 18^>8 
S2 


Then,  using  the  invariant  imbedding  principle,  we  have 


(1)  U^^(x)  -  F^(A)  jl  +  U^(x  -  A)j 


♦  (1  -  *(1)(A)) 


U^(x-A)+4 


Z 

i-4 


1  +  F 


^(A)U^(x  -  A)j 


♦  w(x,A), 

where  Probjw(x,A)  /  oj  -  o(a). 

Let  us  interpret  this  equation.  The  superscripts  in  them¬ 
selves  have  no  significance  and  serve  merely  to  distinguish  one 
random  variable  from  another.  Notice  that  if  the  initiating 
particle  makes  a  collision  in  passing  from  x  to  x  —  A 
(Fig.  13),  the  second  term  of  (l)  is  zero  while  the  first  gives 
Just  the  number  of  particles  emergent  from  x  -  A  — namely,  the 

one  Immediately  out  due  to  the  collision  in  (x  —  A,x)  plus 

(2) 

the  random  number  U'  '  due  to  the  left  moving  neutron  acting 
as  a  source  particle  for  the  rod  (0,x  —A).  If  there  is  no 
Initial  interaction  in  the  Interval  (x  —  A,x),  then  the  first 
term  of  (l)  is  zero,  and  the  second  term  counts.  Of  the  random 
number  (x  —  A)  of  neutrons  emerging  from  x  —  A  some 

make  no  collisions  in  (x  —  A,x)  and  hence  contribute  only  one 
to  the  sura.  Others  make  a  collision  in  (x  -  A,x)  giving  not 
only  an  immediate  right  traveling  particle,  but  also  a  random 
number  U^(x  —  A)  reflected  from  (0,x  —  A).  The  fact  that 
all  other  processes  are  "unlikely"  is  included  in  the  term  w. 
We  now  Introduce  the  generating  function 


(2) 
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u(x,B  )  -  sis U(X)/  . 

i  ; 

Then,  using  the  properties  of  the  generating  function  and 
writing  Prob|F^(&)  -  l>  -  oA  4  o(A),  we  find,  after  care¬ 
ful  calculation  with  (1),  the  equation 

(3)  u(x,a)  ■  dA su(x  -  A, a) 

4  (1  -  dA)u(x  -  A,  (1  -  oA)s  4  0Aau(x  -  A,») )  +  o(A) 

■  cAsu(x,b)  +  (1  -  oA)[u(x  -  A,s) 

+  u  (x  -  A,s)(-  dAs  4  dAsu£]  4  o(A). 
s 

This  leads  at  once  to  the  partial  differential  equation 

(4)  u  «  -  du  +  dsu  4  dsu  4  dsuu  . 

This  is  j  elsely  Equation  (25-5)  derived  there  by  using  quite 
different  methods,  and  from  it  may  be  obtained  the  usual  flux 
equations  for  the  rod  case .  It  was  first  pointed  out  to  us  by 
T.  E.  Harris  that  the  functional  equation  approach  could  be 
applied  directly  to  derivation  of  the  generating  function. 

This  is  an  application  of  a  quite  general  principle  that 
methods  suitable  for  the  derivation  of  first  moments  can  be 
used  almost  unchanged  to  derive  generating  functions .  Prom 
these,  relations  for  the  higher  moments  can  be  readily  obtained. 

See  ^51  for  another  illustration. 

28.  A  Basic  Stochastic  Functional  Equation 

We  shall  now  derive  a  basic  stochastic  functional  equation 
applicable  to  the  generalized  situations  described  in  §26.  We 
introduce  the  appropriate  random  variables . 


Let 
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(l)  11(2,3’, p)  m  random  number  of  particles  In  state  S' 

reflected  from  T)  over  all  time  due  to  an 
initial  particle  in  state  S  impinging  on 
at  time  zero.  (S  and  S’  may  now 
include  information  as  to  the  specific 
location  on  7p. ) 


(2)  r(S,7?) 


1  if  the  particle  in  state  S  is  involved  in 
<  an  interaction  in  the  stratum  {pp  -  A, 7?), 

0  otherwise. 


fi 


(?)  Z(3;S1',S2',...,3k';7?)  ml 


if  the  result  of  an  inter¬ 
action  of  a  particle  in  state 
S  with  the  medium  is  to  pro¬ 
duce  k  particles  in  state 


0 


S} ’ ,Sg ’ ,  •  •  • 
otherwise . 


k 


1,2,5,.. 


By  the  stochastic  variables  ,r^  ,Z^  ^  we  shall  mean 

respectively  any  of  a  denumerable  set  of  variables  with  the 
properties  described  above.  Different  superscripts  will  serve 
merely  to  distinguish  different  random  variables,  and  sometimes 
they  will  be  omitted. 

Alao,  let  Q(S,7?)  denote  the  deterministic  change  in 
state  caused  by  passage  through  the  stratum  (7?  -  A,p) .  By  a 
deterministic  change  we  mean  one  that  occurs  by  virtue  of  mere 
positive  change  without  any  interaction  necessarily  talcing 
place.  Por  example,  any  dependence  of  on  r  mentioned  in 
6 1 7  results  in  such  a  change. 
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Finally,  we  assume  that  only  a  1'lnite  set  of  states  S 
is  possible.  This  permits  us  to  use  characteristic  functions 
rather  than  characteristic  functionals,  and  considerably 
simplifies  the  discussion. 

At  the  risk  of  being  somewhat  redundant,  we  now  discuss 
the  random  physical  processes  which  occur.  A  particle  in 
state  S  Incident  on  the  surface  7?  enters  the  stratum 
(?7  —  A, >7).  There  it  may  undergo  a  deterministic  transformation 
converting  it  from  state  S  to  state  In  addition,  it 

may  undergo  a  stochastic  change  (collision)  producing  a  random 
number  of  particles  in  a  random  set  of  states.  Each  of  these 
particles  acts  independently  of  the  others  and  may  emerge  from 
?7  or  be  incident  upon  7?  -  A.  In  the  latter  case  the  particle 
is  a  "source"  particle  on  77  -  A,  resulting  in  a  random  number 
of  particles  in  a  random  set  of  states  eventually  emerging  from 
77  —  A.  Each  of  these  reflected  particles  may  undergo  deter¬ 
ministic  and  stochastic  transformations  in  (^7— A, ^7).  Some 
emerge  from  >7  and  are  counted,  others  return  as  source  particles 
into  77  —  A  and  must  be  followed  out  again.  No  processes  need 
be  traced  to  order  higher  than  A  — that  is,  to  more  than  one 
collision  in  (7?—  A,  >7). 

Enumerating  these  events  mathematically,  we  arrive  at  the 


functional  equation 
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(4)  U(S,S',?) 

-  r(S,?7)  Z  Z(SjS  •,S?',...,S  ')?) 

[Vj 

+  (1  -  r(S,?l))  Z  (1  -  r(l>(S',?7)) 


z  tr(1)(3.  '.S', y-  a) 
i-1  1 


i-3 


+  (1  -  r(S,?))  Z 

S,  ' 


m  u 


n' 

Z 

i-1 


Hsi  ] 


Z  u(p)(s  ",3',9-  6)j 

P-1  ^ 


+  w(3,S',^, 


where 


n  -  U(Q(S,??),S’,77-  A), 

n'  -  D(Q(S#7?),Sm',^-  A), 

S'  Is  auoh  that  S’  -  Q(3>',??), 

and  w(S,S',??)  la  the  contribution  from  events  that  have 
probability  o ( A) •  The  symbol  js^j  means  all  the  subsets  of 
the  set  of  states. 

Proceeding  as  in  §27,  it  is  now  theoretically  possible  to 
derive  the  flux  equation  of  §27  by  taking  expected  values  of 
(4).  it  is  likely  that  the  direct  approach  of  §26  is  really 
simpler.  However,  the  basic  stochastic  equation  (4)  seems  of 
considerable  theoretical  Interest  and  could  well  be  used  for 
computational  purposes  in  place  of  a  direct  Monte  Carlo  approach. 
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29.  Collision  Processes 

So  far  we  have  not  considered  transport  processes  in  which 
the  particles  interacted  with  each  other,  nor  processes  in 
which  the  medium  was  affected  by  the  transport  process  .  In  the 
sections  that  follow,  we  shall  consider  processes  of  this  nature. 

Let  us  discuss  a  transport  process  in  a  one-dimensional 
rod  in  which  we  allow  neutron-neutron  interactions,  the  result 
being  annihilation.  To  simplify  matters,  let  us  suppose  that 
there  is  no  energy  dependence,  and  that  collisions  occur  only 
between  neutrons  travelling  in  opposite  directions  along  the 
line.  As  usual,  we  shall  suppose  that  when  fission  occurs  one 
neutron  is  produced  in  the  forward  direction  and  one  in  the 
backward  direction. 

Let 

(1)  (a)  6b  +  o(A)  =  the  probability  that  a  neutron  will 

interact  with  a  segment  of  length  A 
and  produce  fission. 

(b)  u(y;x,z)  «*  the  expected  number  of  neutrons  per  unit 

time  passing  an  interior  point  y  to  the 
right,  as  a  result  of  z  neutrons 
per  unit  time  introduced  at  x  (see 
Pig.  14  below). 

(c)  v(y;x,z)  »  the  expected  number  of  neutrons  per  unit 

time  passing  an  interior  point  y  to  the 


left. 
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(d)  x(u,v)A  +  o(A)  ■  the  expected  number  of  neutrons  in 

stream  of  strength  u  which  are 
annihilated  per  unit  time  due  to 
collisions  with  an  opposing  stream 
of  strength  v,  in  an  interval  of 
length  A. 


y-A 


y+A 


Pig.  14 


JfrO .  Intern  1  Plux  Equations 

The  usual  "input— output  **  analysis  yields  the  relations 

(1)  u(y)  -  u(y-A)(l-dA)  +  u(y-A)dA+  v(y)oA 

-  Ak(u,v)  +  o(A) , 

v  (y )  *  v  (y  +  A)U  -  +  v(y)oA  +  u(y)dA 

-  Ak ( u , v )  +  0  ( A) , 

where  we,  at  this  time,  suppress  the  dependence  upon  x  and 
z,  and  write 

(2)  u(y,x,z)  b  u(y),  v(y,x,z)  ■  v(y). 

Passing  to  the  limit  as  A  — *  0,  we  are  led  to  the  nonlinear 
system  of  differential  equations 


(3) 


u* (y)  ■  dv  -  k(u, v) , 
v *  ( y )  «  -  du  +  k(u, v) . 
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The  boundary  conditions  are 
(4)  u(0 )  -  0,  v{x)  -  z, 

two-point  boundary  conditions  . 

51.  Discussion 

In  general,  the  equations  in  (>0.5)  cannot  be  resolved  in 
terms  of  elementary  transcendents  of  analysis.  Consequently, 
if  we  wish  to  obtain  a  numerical  solution  of  (>0.>)  cum  (>0.4), 
we  must  resort  to  various  computational  schemes.  Although  a 
number  of  these  are  available,  it  cannot  be  said  there  are  any 
of  automatic  application.  Questions  of  this  nature  are  of 
great  difficulty,  and  perhaps  are  most  easily  handled  by  being 
bypassed.  In  the  next  section  we  shall  approach  this  problem 
in  a  different  way. 

If  we  make  the  assumption  that  k(u,v)  ■  buv,  b  >  0,  a 
certain  amount  of  analysis  can  be  carried  out.  See  [>4]  for 
some  theoretical  and  computational  results.  Henceforth,  we 
assume  this  form  of  k(u,v). 

>r. .  Reflected  and  Transmitter  Flux 

Let  us  now  approach  tKis  problem  by  means  of  functional 
equation  techniques.  Let 
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(l)  r(z,x)  -  the  expected  number  of  neutrons  reflected 

per  unit  time  from  a  homogeneous  bar  of 
length  x  as  a  result  of  having  z  neu¬ 
trons  incident  at  x  per  unit  time. 


♦  r(z,x) 


0  x 

Fig.  15 

To  evaluate  the  expected  number  of  neutrons  reflected 
from  a  bar  of  length  x  +  A  we  note,  first  of  all,  that  some 
collisions  with  nucleil  may  occur  immediately  when  the  z 
neutrons  enter  the  segment  [x,x  +  A]  .  Since  each  such  col¬ 
lision  results  in  a  neutron  going  to  the  right,  dzA  particles 
emerge  at  x  +  A.  Meanwhile,  since  a  neutron  is  produced  going 
to  the  left  also,  the  original  flux  z  is  not  affected  by 
collisions  with  nucleil.  However,  this  flux  is  reduced  by 
annihilation  by  an  amount  bzr(z;x)A  due  to  the  opposing  flux 
out  of  x.  It  is  also  Increased  by  an  amount  dr(z;x)A  due  to 
fission  collisions  in  [x,x  +  A]  made  by  the  flux  out  of  x. 

Hence  there  is  a';  x  a  source  of  strength  z  -  bzr(z;x)A  + 
d r(z;x)&.  Finally,  the  reflected  flux  due  to  this  source  is 
partially  annihilated  by  interactions  with  the  impinging  flux 
in  [x,x  +  A]  .  Summing  up,  we  are  led  to  the  relation 

(2)  r(z;x  +  A)  »  dzA  +  r(z  —  brzA  +  dAr(z,x) )  [l  -  bzA]  ♦  o(A)  . 
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By  letting  k  tend  to  zero  vfe  find  that  r(z;x)  satisfies  the 
quasilinear  first  order  partial  differential  equation 

(3)  r  «  dz  -  bzrr  +  drr  -  bzr 

where,  as  usual,  the  subscripts  indicate  partial  differentiation. 
The  reflection  function  r(z;x)  also  satisfies  the  initial 
condition 


(4)  r(z;0)  -  0  and  r(0;x)  -  0. 


The  Equation  (3)  specializes,  for  b  «  0,  to  the  Riccati 
equation  derived  in  earlier  sections  for  the  reflection  coef¬ 
ficient  (§4).  It  may  be  resolved  via  characteristic  theory 
[55]  or*  by  direct  numerical  integration,  returning  essentially 
to  (2).  The  equations  for  the  characteristics  are 


(5) 


dx 

3s  “ 


1, 


as  '  bzr  -  dr> 

=  dz  -  bzr. 


Since  x  =*  s,  z-0,  r  *  0  is  a  solution  of  the  system  (5) 
passing  through  the  point  x  -  z  -  r  ■  0,  we  find  that 


(6)  r(0;x)  -  0, 


as  was  assumed  above  on  physical  grounds. 

Once  the  function  r(z;x)  has  been  determined  for  suitable 
ranges  of  z  and  x,  one  may  reduce  the  determination  of  u(y) 
and  v(y),  the  Internal  fluxes,  to  the  solution  of  initial  value 
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problems,  as  was  mentioned  earlier.  If  the  incident  flux 

v(x)  ■  z  is  specified,  then  the  reflected  flux  is  r(z;x)  -  u(x), 

so  that  now  both  u(y)  and  v(y)  are  specified  at  y  »  x. 

Through  use  of  Equations  (30.3)  the  functions  u(y)  and  v(y) 
may  now  be  determined  on  the  entire  interval  [o,x]  . 

The  equations  satisfied  by  the  transmitted  flux  t(zjx), 

where 

(7)  t(z;x)  -  the  expected  number  of  neutrons  emergent 

from  the  end  y  -  0  of  a  homogeneous  bar 
of  length  x  as  a  result  of  having  z 
neutrons  per  unit  time  incident  on  the  end 
y  -  x, 

are  similarly  derived.  We  have 

(8)  tx  -  (<J  -  bz)rtz 

along  with  the  boundary  conditions 

(9)  t(Ojx)  -  0,  t(z;0)  -  z. 

33.  Discussion 

Numerical  solution  of  the  foregoing  Equation  (32.3) ,  can 
be  obtained  either  by  use  of  (32.2),  by  means  of  conventional 
techniques,  or  by  means  of  the  characteristics. 

In  [34]  will  be  found  a  brief  discussion  of  collision 
processes  of  this  nature  with  energy  dependence.  In  the 
remainder  of  this  paper  we  consider  only  processes  without 
particle— particle  interaction. 


34.  Time  Dependent  Rod  Case — Internal  Flux 
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Ihus  far  we  have  concentrated  our  attention  on  the 
stationary  state  in  dealing  with  both  internal  and  reflected 
and  transmitted  fluxes.  In  this  section  we  shall  obtain  the 
equations  for  the  internal  flux  in  the  rod,  taking  time  vari¬ 
ation  into  account.  We  assume  that  the  rod  has  all  the  usual 
properties,  and  that  neutrons  in  It  travel  with  speed  c. 
Thus,  a  time  A/c  is  spent  in  traversing  a  distance  A.  It 
is  now  relatively  easy  to  see  tt  at  if 

(1)  u(y,t)  -  average  number  of  neutrons  passing  y  per 

second  at  time  t  and  going  to  the  right; 
v (y,t )  -  average  number  of  neutrons  passing  y  per 
second  at  time  t  and  going  to  the  left; 


then 


(2)  u (y  +  A,t  +  |)  *  u(y,t)(l  -  dA)  +  dAu(y,t) 

+  v(y  +  A, t )oA  +  o(&), 

v(y,t  +  | )  -  v (y  +  A,t)(l  -  dA)  +  dAv(y,t) 
+  u(y ,  t  )dA  +  o  (A) , 


giving  in  the  limit  as  A 


(3) 


1  9u 

c  £T  5y 


cfv, 


1  2v  7w  . 
c  ~Sl  -  Jy  “  du- 


0  the  equations 


Equations  (3)  are  subject  to  the  boundary  and  initial 


conditions 
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(M  u(0,t )  -  0, 

v(x,t)  -  1, 
u(y,0)  -  v (y ,0)  -  0, 

for  the  case  of  the  rod  with  a  unit  flux  at  the  right  end, 
y  -  x,  imposed  at  t  =  0. 

This  formulation  is  the  classical  one.  Equation  (If.*) 
already  mentioned  without  discussion  is  the  analogue  of  (3) 
for  the  more  complicated  geometry.  Such  equations  as  (15*1) 
are  readily  derived  using  the  principles  of  this  section. 

35  •  Time  Dependent  Rod  Case — Reflected  Flux 

.now  turn  to  a  formulation  of  the  time  dependent  rod 
problem  using  the  basic  ideas  of  invariant  imbedding.  Let  us 
consider  a  single  "trigger"  neutron  entering  the  rod  at  x 
at  time  t  «  0.  It  is  uhen  convenient  to  Introduce 

(1)  U(x,t)  ■  total  number  of  neutrons  reflected  from  x 

up  to  time  t  due  to  the  onef  trigger 
neutron  in  at  x  at  time  t  -  0. 

Clearly 

(2)  U(x,t)  «  P%  u(x,s)ds, 

% 

where  u  is  the  function  of  §34,  now  subject  to  the  delta 
function  type  initial  condition.  We  propose  to  imbed  the  rod 
of  length  x  as  usual  (see  Figure  16). 
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Pig.  16 

It  should  be  noted  that  this  also  provides  an  Imbedding  in 
time,  and  that  the  important  time  increment  is  2A/c,  repre¬ 
senting  the  time  it  takes  for  a  particle  to  cross  [x,x  +  A] 
in  each  direction. 

We  find  as  usual  that  there  may  be  an  immediate  collision 
in  [x,x  +  A]  ,  providing  a  flux  out  to  the  right  of  dA. 

There  is  also  a  single  trigger  neutron  into  [0,x] ,  and  this 
enters  at  time  t  *  A/c .  At  time  s  +  A/c  let  there  be 
u(x,a)  neutrons  per  second  emergent  at  x.  Clearly 

(5)  u(x,s)  -  (x,s) . 

Some  of  these  neutrons  make  collisions  in  [x,x  +  A ] . 

There  is  surely  a  flux  out  of  x  +  A  of  rate  u(x,s).  How¬ 
ever,  half  of  the  fission  neutrons  return  to  x  at  this  time. 

These  provide  a  new  source  into  x,  and  will  contribute  to 
the  total  flux  out  during  the  remaining  t  -  s  units  of 
time.  Thus 

(4)  U(x,t  +  —■)  «  6h  +  U(x,t)  +  dA  u(x,s)U(x,t  -  s)ds  +  o(A) 

c 

giving  the  equation  of  mixed  differential— integral  type 

+  I  ft  "  0  +  u(^»s)U(x,t  -  s )ds , 


(5) 
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subject  to  the  conditions 

(6)  U(x,0)  -  0,  U(0,t)  -  0. 

Wie  system  (5)  —  (6)  has  been  analyzed  rigorously  (see 
[45]  ) .  Ihe  convolution  form  of  the  integral  tern  makes 
possible  an  explicit  analytic  representation  of  the  solution 
when  the  Laplace  transform  is  used. 

56.  Modification  of  Medium  during  Transport  Process 

In  all  that  has  gone  before,  we  have  supposed  that  the 
properties  of  the  medium  have  remained  unaltered  by  the  trans¬ 
port  process  occurring  within  it.  This  Is,  of  course,  always 
an  approximation  of  greater  or  lesser  validity. 

Two  very  interesting  processes  in  which  interaction  with 
the  medium  are  taken  into  account  are  the  free  boundary  problems 
of  hydrodynamics,  and  the  Stefan  problems  of  heat  conduction. 

As  a  preliminary  to  an  investigation  of  Stefan— type  processes 
by  invariant  imbedding  techniques,  we  wish  to  discuss  a  one- 
dimensional  transport  process  in  a  rod  whose  length  is  changing 
as  a  function  of  time. 

57.  The  Physical  Process  and  its  Mathematical  Formulation 
Let  us  consider  a  rod  which  extends  from  0  to  x  at 

time  t  -  0.  The  rod  grows,  or  erodes,  at  a  specified  rate  so 
that  th^  position  of  the  left  end  is  given  by  *  f(t), 

f(0)  »  0,  while  the  right  end  remains  fixed,  ■  x.  (See 

the  figure  below.)  A  neutron  traversing  a  distance  A  in  the 
rod  has,  as  before,  probability  dA  +  o(A)  of  suffering  a 
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collision  with  elements  of  the  medium.  In  the  event  of  a 
collision,  which  produces  fission,  two  neutrons  emerge,  one 
moving  to  the  left,  the  other  to  the  right.  No  neutron  can 
re-enter  the  rod  once  having  left  it,  and  all  neutrons  have 
constant  velocity  c.  A  single  neutron,  the  "trigger,"  enters 
the  rod  at  x  at  time  t  -  s .  We  ask  for  the  expected  total 
number  of  neutrons  that  emerge  at  the  right  end  up  to  time  t, 
denoting  this  quantity  by  U(x,s,t). 

It  must  be  noted  that  the  condition  (—  f'(t))  <  c  for 
all  t  is  convenient  to  prevent  neutrons  from  being  rather 
artificially  "trapped"  in  the  rod.  We  shall  not  discuss 
interesting  problems  of  this  nature  here  . 

To  apply  the  invariant  imbedding  method,  we  immerse  the 
original  process  within  the  class  of  all  processes  indexed  by 
x  >  0,  and  then  express  the  relationship  between  neighboring 
processes.  The  left  ends  are  uniformly  to  obey  the  same  law, 

-  f(t),  as  that  of  the  original  rod. 


*L-° 

t  =  0 

_  j — - -  - — . . 1  « -  Neutron 

XL-f(§  -  a/c)  XR=x  x+A 

t  «  s  -  A/c 


XL-f(t)  0  XR-x  x+A 

t  >  • 

Pig.  17.  Neutron  transport  in  a  rod  of  varying  length. 


v* 


x+A 
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Let  us  now  analyze  the  process.  At  time  t  ■  s  —  A/c  the 

trigger  neutron  enters  at  x  4-  A,  and  may  or  may  not  suffer  a 

collision  in  passing  from  x  +A  to  x.  If  it  does,  a  single 

neutron  immediately  emerges  at  the  right.  In  either  event,  a 

neutron  passes  xR  at  time  t  -  s.  This  acts  as  a  trigger  for 

the  original  rod,  and  produces  a  flux  emergent  from  x^  at 

times  t  >  s.  The  expected  rate  of  emergence  of  such  neutrons 

on 

is  given  by  u(x,s,t)  «  ^-(x ,s , t )  .  Neutrons  in  this  flux  may 
or  may  not  suffer  collisions  with  the  rod  material  in  going 
from  x  to  x  +  A.  In  either  event  U(x,s,t)  neutrons  emerge 
from  x  4-  A  by  time  t  4-  A/c  .  In  addition,  those  neutrons 
which  do  make  collisions  in  (x,x  4-  A)  at  time  t',  s  <  t'  < 
t  4-  A/c,  contribute  trigger  neutrons  at  at  a  rate 

dAu(x,s,t')  4-  o(A).  The  resultant  flux  out  of  x  4  A  up  to 
time  t  4-  A/c  is  then  dA  u(x,3,t,)U(x,t,,t)dt'  4-  o(A). 

o 

All  other  processes  yield  contributions  of  order  o(A).  Thus 
we  find 

(1)  U(x  4-  A, 8  -  £,t  4  |)  -  dA  4-  U(x,s,t) 

4  dA^/>tu(x,s,t'  )U(x,t '  ,t )dt’  +  o(A). 

8 

Hence,  letting  A  — ►  0,  we  obtain  the  equation 

(2)  +  If?  •  *  +  <5/t  u(x,a,t')U(x,t',t)df. 

The  conditions  imposed  on  U  are 
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(3)  U(x,8,t)  -  0,  »  2 

U(x,s,t)  -  0,  f(s)  -  XR. 

•Rie  first  merely  states  that  no  neutrons  emerge  by  time  t  if 
the  trigger  enters  at  a  later  time,  while  the  second  reflects 
the  fact  that  none  can  emerge  to  the  right  if  the  rod  is  of 
length  zero  when  the  trigger  neutron  impinges. 

38.  Discussion 

Equation  (37*2)  has  the  same  general  structure  as  that 
studied  in  our  earlier  sections  for  the  case  of  the  rod  of 
fixed  length  with  the  neutron  entering  at  t  -  0.  It  is  not 
difficult  to  reduce  (37*2)  to  the  simpler  equation  when 
f(t)  =  0 .  However,  the  fact  that  the  Integral  term  in  (.37.2) 
is  no  longer  of  conv  ^ution  form  Introduces  new  complications 
in  studying  propei  cs  of  the  solution. 
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Part  III 

DIPFUSION  THEORY — A  LIMITING  CASK 
OP  TRANSPORT  THEORY 

39.  Diffusion  aa  a  Limiting  Process 

In  previous  sections  in  this  paper,  we  have  investigated 
a  variety  of  simple  models  of  transport  theory  by  means  of  the 
functional  equation  technique  of  invariant  imbedding. 

Neutrons  are  mathematically  abstracted  to  be  point  particles 
with  finite  velocities,  while  fission  and  scattering  are 
characterized  by  certain  probabilities  (cross  sections)  of 
branching  and  reversal  or  reorientation  of  direction  in  the 
medium  within  which  the  process  is  occurring.  In  the  great 
proportion  of  cases  we  assume  no  neutron-neutron  interaction, 
and  no  change  in  the  properties  of  the  medium  over  time,  al¬ 
though  we  have  discussed  both  of  these  phenomena  to  slight 
extents,  ^29  and  ^37. 

It  la  of  interest  for  several  reasons,  from  both  the 
mathematical  and  physical  points  of  view,  to  discuss  in  detail 
what  happens  to  the  various  categories  of  transport  equations 
derived  from  different  applications  of  invariant  Imbedding  as 
the  velocity  of  the  neutron  is  allowed  to  become  arbitrarily 
large  with  a  corresponding  increase  in  the  probability  of  a 
collision . 

This  idea  is  a  quite  natural  one  and  one  that  has  been 
pursued  by  a  number  of  different  investigators  with  di  tb^ent 
alms  in  mind.  Diffusion  theory  classically  has  bee  re*,  w-ued 
as  an  approximation  to  the  more  rigorous  (but,  ol  v^urse  not 
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completely  rigorous)  transport  theory  under  the  assumption  of 
high  velocity  and  small  mean  free  path  [4o]  .  Furthermore, 
passage  to  the  limit  in  the  "telegrapher's  equation,"  a  linear 
partial  differential  equation  of  hyperbolic  type,  has  been 
carried  out . 

Prom  another  direction,  the  discrete  random  walk  process 
yields  the  diffusion  equation  in  the  limit.  This  observation 
has  been  made  the  basis  for  a  considerable  amount  of  analytic 
and  computational  effort,  centering  about  the  theme  of  "Monte 
Carlo"  techniques. 

Our  principal  aim  here  is  to  study  the  limits  of  the  non¬ 
linear  functional  equations  obtained  from  the  transport  pro¬ 
cesses  with  finite  velocity,  e.g.  those  appearing  in  §54  and 
§55,  as  the  velocity  increases  without  bound.  In  this  way,  we 
obtain  corresponding  results  for  heat  or  diffusion  processes, 
where  the  physical  picture  is  not  as  clear.  Having  obtained 
the  equations  in  this  indirect  and  complex  faahlan,  we  can 
then  interpret  them  in  such  a  way  as  to  be  able  to  derive  them 
directly  by  invariant  imbedding  techniques.  In  all  cases,  the 
equations  are  of  the  generalized  Rlccatl  type  which  we  recog¬ 
nize  as  characteristic  of  these  processes  of  mathematical 
physics,  when  described  in  invariant  imbedding. 

At  the  present  time,  we  are  studying  the  question  of 
treating  Stefan— type  diffusion  problems  by  a  similar  passage 
to  the  limit  in  the  equations  derived  from  transport  processes 
with  variable  boundaries.  This  is,  as  might  be  expected,  a 
complex  problem.  Some  initial  results  have  been  given  above,  §57* 
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Throughout  the  sections  that  follow,  we  shall  use  a  simple 
generalization  of  the  idealized  one-dimensional  rod  process 
treated  in  the  foregoing  pages.  In  the  following  section,  we 
shall  obtain  some  new  equations  for  the  flux  within  the  rod, 
assuming  finite  velocities  initially.  In  §4l,  we  derive  Pick's 
law  for  this  simple  process .  This  is  important  for  our  pur¬ 
poses,  since  It  is  the  analysis  of  this  result  which  suggests 
the  combinations  of  functions  which  should  be  used  in  the 
limiting  case.  In  §42,  we  study  the  limiting  form  of  the 
internal  flux  as  the  velocity  becomes  infinite,  and  in  §4}  the 
diffusion  process  giving  rise  to  the  function  obtained  in  this 
way  is  analyzed. 

We  then  turn  to  our  primary  objective,  the  passage  to  the 
limit  of  the  nonlinear  integro— differential  equation  obtained 
for  the  reflected  flux  in  the  neutron  transport  transport  case 
by  means  of  the  technique  of  invariant  imbedding.  In  §45,  we 
show  how  to  obtain  the  result  by  direct  application  of  the 
imbedding  technique  to  the  diffusion  process. 

Throughout  this  part  of  the  paper,  our  methods  are  again 
largely  formal,  since  we  are  principally  interested  in  demon¬ 
strating  J.e  applicability  of  invariance  principles  .  The 
exis4  ce  of  relevant  limits  and  the  applicability  of  Laplace 
transform  methods  are  taken  for  granted  in  order  to  arrive 
iuiOH-y  at  the  desired  equations.  These  questions  can  be 
tvadied  in  a  rigorous  fashion,  and  for  the  simple  mathematical 
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models  considered  here,  there  is  little  difficulty  in  carrying 
out  this  program.  Since,  however,  we  know  that  the  passage  to 
the  limit  involves  a  reduction  from  a  hyperbolic  partial  dif¬ 
ferential  equation  to  a  parabolic  partial  differential  equation, 
Involving  inter  alia  a  redundancy  in  the  initial  conditions,  we 
can  expect  some  difficulties  in  the  general  case.  The  corres¬ 
ponding  study  for  ordinary  differential  equations  when  a 
limiting  value  of  a  parameter  results  in  a  drastic  change  in 
the  order  of  the  equation  is  of  some  subtlety.  Particularly 
Interesting  examples  of  equations  of  this  nature  occur  in 
various  hydrodynamical  investigations  where  viscosity  plays 
the  role  of  the  parameter  which  approaches  zero;  of.  Wasow  [4l]  . 

40.  The  Transport  Equation 

To  begin  our  work  it  is  necessary  to  write  down  some  trans¬ 
port  equations  in  fairly  general  form.  While  some  of  them  are 
not  to  be  found  in  the  literature,  they  may  be  readily  derived 
by  the  methods  of  previous  sections. 

Consider  a  rod  of  material  which  transports  neutrons,  and 
let  the  neutrons  have  constant  velocity  j  (raonoenergetic 
case).  The  usual  collision  processes  take  place  with  the 
probability  of  a  collision  in  a  length  A  of  the  rod  taken 
to  be  oA  +  o(&),  where  6  is  a  constant.  On  the  average, 

2k  neutrons  emerge  from  a  collision  inside  the  rod,  k  going 
to  the  left  and  k  going  to  the  right.  We  take  the  rod  to 
extend  from  0  to  x  (see  Pig.  18),  and  designate  the 
coordinate  of  an  interior  point  by  y. 
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I - 1 - H - 

0  y 

Pig.  18 

To  initiate  the  process,  we  suppose  that  there  is  a  tine 
dependent  source,  q(t)  neutrons  incident  to  the  left  per 
second  at  x,  and  none  at  the  end  0.  Finally,  we  suppose 
that  particles  emergent  at  0  or  x  cannot  re-enter  the  rod. 
We  write 

(1)  u(y,t)  =  the  average  number  of  neutrons  per  second 

passing  y  at  time  t  and  moving  to  the 

right, 

v(y,t)  *  the  average  number  of  neutrons  per  second 
passing  y  at  time  t  and  moving  to  the 

left . 

Using  the  methods  outlined  in  ^34,  it  is  easily  found  that 

(2)  fy  f  a  y  dkv' 

~  %  +  \  %  *  tfku  +  6{k  “  1)v' 

u(0,t)  »  0,  v ( x , t )  =  q(t),  t  >  0; 

u  { y ,  0 '/  »  v  (y  ,0)  -  0,  0  <  y  <x. 

For  some  purposes  it  is  convenient  to  talk  about  the  total 
flux  from  time  0  to  t.  We  shall  consistently  use  capital 
letters  to  Indicate  quantities  integrated  over  time.  Thus 


-j  4 —  Source  q(t) 

x 
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(3)  U(y,t)  -  /)t  u(y,z)dz, 

‘/0 


Q(t)  -  /}t  q(z)dz, 
^0 


etc.  It  is  easy  to  see  that  the  integrated  quantities  satisfy 
equations  identical  to  (2)  with  the  lower  case  letters  being 
replaced  by  capitals. 

At  times  it  will  be  desirable  to  make  clear  the  type  of 
source  in  the  particular  problem  under  discussion.  Hence  we 
shall  occasionally  write 

(4)  u(y,t;q), 

U(y,t;q) 

etc.,  to  emphasize  the  source  in  question.  The  source  will  be 
deleted  when  the  meaning  is  clear.  In  these  cases  we  shall 
write  u(y,t;q)  •  u(y,t),  and  so  on.  In  other  instances  the 
source  may  be  indicated  and  the  dependence  upon  y  or  t 
left  out. 

Consider  now  the  case  in  which  the  source  consists  of  a 
single  "trigger"  neutron  at  t  »  0.  Thus,  formally, 
q ( t )  «  &(t),  where  b  is  the  Dirac  delta  function.  We  focus 
attention  on  the  particles  reflected  from  the  rod  at  x, 
writing  r(x,t;8)  for  the  number  emergent  per  second,  and 
R(x,t,S)  for  the  total  number  oraergent  up  to  time  t. 

Clearly 


(5) 


r(x,t;b)  =  u( x, t ; b  ) , 
R(x, t ; 8  )  =  U(x,t;b  ) . 
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However,  it  Is  again  well  to  regard  the  x  in  the  arguments 
of  r  and  R  as  referring  to  the  length  of  the  rod  rather 
than  to  the  coordinates  of  the  end  point  of  the  rod.  With 
this  rather  subtle  distinction  in  mind  one  then  finds,  pro¬ 
ceeding  as  oefore  (^35), 

(c)  +  ok  )R(fe) 

+  ak  /?t  r(x,z;6)R(x,t  -  z;8)dz, 

60 

R(x,0;O  )  ---  R(0 ,  t ; 8  )  =  0. 

Notice  that  this  characterizes  the  reflected  flux  in  a  fashion 
Independent  o f  the  Internal  fluxes . 

For  the  corresponding  case  in  which  there  is  a  source 
q(t),  we  note  that  the  fundamental  physical  process  is  addi¬ 
tive,  as  a  consequence  of  our  tacit  assumption  that  there  are 
no  interactions  between  neutrons  passing  in  opposite  directions, 
we  can  write 

(7)  R ( x , t ; q )  -  /?t  q(z)R(x,t  -  z;8)dz. 

°0 

To  find  an  equation  satisfied  by  R(x,t;q),  we  utilize 
the  bap  lace  transform,  writing 

R(x,s)  e-st 

60 


(&) 


R(x, t )dt , 
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with  a  consistently  similar  notation  for  transforms  of  other 
functions.  Then,  from  (6),  with  q  =  b, 

dRr(b)  o  O 

(9)  — jfc—  +  §  sRl(&)  =  x  +  2a(k  -  l  )rl(6)  +  oksRL2(&), 

Rr  (0,8;S)  =  0. 

ij 

From  (7), 

(10)  Rjx,s;q)  =  qL(s)RL(x,s;fe). 

Hence , 

dRr  (  Q  )  Q 

(11)  +  -  3RL(q)  =  T  qL  +  20(R  -  l)RL(q) 

+  aksRL(q)RL(&) , 

which  leads  back  to 

(12)  +  |  =  okQ(t)  +  2o( k  -  l)R(q) 

+  ok  p 1  r(x,z;6)R(x,t  -  z;q)dz, 

c0 

R(0, t;q)  =  R(x, 0; q )  =  0. 

This  clearly  reduces  to  (6)  when  q(t)  =  &(t). 

We  shalL  derive  one  other  special  case  of  (12),  corres¬ 
ponding  to  the  case  when  q(t)  =  1.  While  conceptually  this 
may  be  a  bit  more  difficult  to  consider  than  the  single  trigger 
neutron  case,  it  has  the  mathematical  advantage  of  avoiding  the 
6- function.  For  this  type  of  source  we  write  the  integrated 
flux  as  R(x,t;l)  or  R( 1 )  and  (12)  becomes 
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(13)  +  f  -  okt  +  2o(k  -  l)R(l) 

+  ok  p  ^  r(x,z;8  )R(x,t  —  z;l)dz. 

^0 


But,  from  (7), 

(14)  R(x,t;l)  =  f)X  R(x,t  —  z;6)dz. 

60 

Then  we  easily  find 

(15)  +  |  =  okt  +  ?r(k  -  i)r(  l ) 

+  ok  pX  r(x,  z;  1  )r(x,  t  —  z;l)dz, 

^0 

R(C,t;l)  =  R(x, 0; 1 )  0. 

4l .  Flck 1  a  Law 

If  we  subtract  the  second  equation  of  (40.2)  from  the 
first  we  obtain 

(1)  -^(U  +  V)  4  ^(U  -  V)  n  -  0(U  -  V). 

For  large  c,  we  expect  the  second  term  on  the  left  to  be 
small.  Hence  we  formally  obtain  the  relation 

(2)  ~{u  +  v)  =  -  o(u  -  v) 

In  the  "limit  of  large  velocity."  Equation  (2)  Is  ordinarily 
referred  to  as  Flck 1 s  Law  [4o]  ,  which  states  that  the  net  flux 


is  proportional  to  the  gradient  of  the  concentration  and  in 
the  opposite  direction. 
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42.  The  Limiting  Case  Obtained  Directly 

To  ootaln  preliminary  results  we  take  Laplace  transforms 
of  (40.2).  Thus,  using  the  notation  introduced  in  (40.7), 

dUT  s 

f1)  37"  +  c  UL  "  “  ^UL  +  <5kvL' 

dvT 

~  +  c  VL  “  dkuL  +  d<k  “  l)vL' 
uL(0,s)  -  0,  vL(z,s)  -  qL(s). 


After  rather  extensive  but  rudimentary  calculations,  we  arrive 
at  the  relations 


(2) 


uL(y,s) 

vL(y,s) 


kCq^(s)  sinh  Xy 

jx  cosh  Xx  +  (2-  +  (1  —  k)o)sinh  Xx} 
qL(a)^X  cosh  Xy  +  .2.  +  (l  -  k)d)sinh  Xy] 

JX  cosh  Xx  +  (2.  +  (l  —  k)d)sinh  Xxj 


where 

(3) 


(1  -  2k). 


tie  now  choose  k  =  1/2,  which  means  physically  that  an 
average  collision  gives  rise  to  one  neutron.  This  choice 
eliminates  the  last  term  in  (3).  Since  we  seek  a  diffusion 
type  equation,  we  let  c  — *  cd  ,  which  is  to  say,  we  allow  the 
velocity  to  become  infinite.  Clearly,  to  preserve  the  process 
we  must  then  require  that  0  — ♦  00  in  such  a  way  that 
lim  c/o  =  D,  a  constant.  Hence,  from  what  has  proceeded, 

lim  ?\  -  Vs/D. 

(it  should  be  noted  that  a  somewhat  more  general  result 
could  have  been  obtained  by  requiring,  instead  of  k  =  1/2, 
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2 

that  11m  o(l-  2k)  =  a.  By  so  doing  we  could  have  accounted 
for  cases  of  absorption  or  fission.  To  do  this  here  would 
merely  complicate  the  ensuing  calculations.) 

Bearing  (50.2)  In  mind,  we  set 


co 


JL(y.»)  -  <J(uL(y,a)  -  vL(y,o) ) , 


J0,L(y's) 


lire 
0-— *oo 


JL(y#s). 


Let  us  consistently  reserve  the  subscript  zero  to  refer  to 
quantities  in  the  limit  as  c  — ♦  00 . 

We  then  discover  that 


(5) 


J0  r(y»s)  “  -  111" 


Oq^(s)X  cosh  Xy 


6 — »oo  j^X  cosh  Xx  +  (-  +  ^)slnh  Xxj 


-  -  2qL(s)  3L- 


s  cosh 


(y^p  lg )  . 


slnh 


(x\/d  1 8  ) 


43.  A  Classical  Diffusion  Problem 

We  now  seek  an  ordinary  diffusion  problem  which  gives  rise 
to  the  limiting  expression  found  in  (41.5).  It  Is  readily  veri¬ 
fied  that  If  Q(y,t)  Is  Implicitly  determined  by  the  relations 


(1) 


2  ^ 

D  =  |l»  0(0. t)  -  0,  0(x,t)  =  2q(t),  9(y,0)  =  0, 


then,  explicitly, 

2q,  ( s  )sinh(>VD-18 ) 

(2)  or(y,a)  -  — ± - - - 

slnh  (xVD  s) 


and 


P-l8^8 


(3)  /L  .  2q.(a)V'lF17  cosh(?^  ■  -»)  , 

dy  ,rzrr 

8inh(xvD  s) 

We  may  summarize  our  results  thus  far  as  follows: 

If  we  consider  the  transport  problem  formulated  in  (40.2) 
in  the  limiting  case  where  c  — *  oo  ,  c/c  — ■>  D,  with  k  =  1/2, 
then  the  problem  is  formally  equivalent  to  the  classical  dif¬ 
fusion  proa lem  ( 1 ) .  The  quantity  lim  (u(y,t)  -+  v(y,t))  may 

o— *oo 

be  Identified  with  ©(y,t),  while  lim  o(u(y,t)  —  v(y,t)) 

. .  .  CH-^CD 

corresponds  to  —  9Q/9y . 

It  Is  possible  to  Identify  ©(y,t)  with  the  total  neutron 
flux  (see  [4o]  )  although  the  dlf  fusion  may  refer  as  well  to 
heat  or  material  concentration.  The  fact  that  a  source  of 
2q(t'>  Is  required  as  part  of  the  initial  conditions  in  the 
problem  (l)  may  De  rather  puzzling  until  one  notes  from  (42.2) 
that,  formally,  both  u(x,t)  and  v(x,t)  approach  q(t)  as 
c  — *  oo  . 

44.  The  Reflected  Flux 

Let  us  now  turn  to  Equation  (40.15)  and  try  to  carry  out 
the  same  type  of  passage  to  the  limit.  It  is  clear  that  we 
must  begin  by  Investigating  the  quantity 

(1)  H(  x,  t ;  q )  -  ajR(x,t;q)  -  Q(t)!, 

'  J 

which  reduces  to  o^R(x, t;l)  —  tj,  when  q(t)  =  1.  Thus, 

(?)  R(l)  =  +  t, 

/  .  s  h( 1 ) 

r(  1 )  =  f  1  • 
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SubBtttuting  these  in  (40.15)  with  k  -  1/2,  we  find 


(3) 


0  /?  t(h(x,  z;  1 )  ,)  (  h(x,t 

2^  1  0  Ml 


—  z 


From  this  we  readily  get 


(4) 


H(x,0; 1 )  =  0,  H(0,t;l)  -  -  at . 


Passing  to  the  limit  as  in  ^43  we  get  (at  least  formally) 

5Hn( 1 )  ,  , 

(5)  ~TT~  +  2D  =  icr  h0(x,z;l)hQ(x,t  -  z;l)dz, 

Hq(x,0; 1 )  =0,  HQ(0,t;l)  =  -00; 

where  HQ(x,t;q)  =  Urn  H(x,t;q),  etc.,  as  agreed.  That  (5) 

0—+00 

is  the  correct  limiting  form  may  be  established  by  a  Laplace 
transform  argument  similar  to  that  of  the  last  section.  We 
omit  the  details. 

It  is  of  some  interest  to  evaluate  HQ.  This  may  be  done 
by  solving  (5),  of  course.  However,  it  is  easier  for  us  to 
note  from  (42.5)  that 

(6)  h0(L(x,s;q)  -  -  gqL(a 

sinh  (xVd  ^s) 

= - —  coth  (x^), 

VsE 


when  q(t)  =1.  We  find 
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(7) 


hQ(x, t ;  1 )  = 


V^rt  D  | 


Y  2  2 

_  xn  \ 

1  +  2  Z  e  r 


n=l 


j 


2  ©  (±  — L_ ) 

x  V2'Dx2j' 


where  0^  la  a  theta  function. 

The  analogue  of  (5)  may  be  derived  easily  for  the  case  In 
which  there  is  an  arbitrary  source  Q(tj.  The  result  is 


(8) 


2Hc(q) 

2x 


+ 


2D  1q(t ) 


z;q)dz, 


H  (x,0;q)  =0,  HQ(0,tjq)  =-00. 


e  now  readily  see  the  following  result: 

If  we  consider  the  transport  problem  formulated  In  (40.2) 

In  the  limiting  case  then  the  quantity  HQ(x,t;q)  la  formally 

equivalent  to  the  quantity  —  ~^~t/  Q(y,z)dz|  where  0 

"y(/C  '  y=x 

Is  defined  by  (43.1).  Further  H^( q )  satisfies  (8)  with 
hQ(l)  given  by  (5). 

45.  A  Direct  Invariant  Imbedding  Approach  In  Diffusion  Theory 
The  equations  thus  far  obtained  are  not  new,  though  our 
approach  to  them  may  be  somewhat  novel.  To  conclude  our  work 
here  we  shall  present  a  method  of  obtaining  (44.8)  by  Invariant 
imbedding  techniques  without  venturing  outside  the  confines  of 
ordinary  diffusion  theory.  The  method  described  holds  promise 
of  being  applicable  In  much  more  complicated  diffusion  processes 
than  that  described  here,  and,  in  particular,  may  eventually 
yield  new  formulations  of  Stefan- type  problems. 
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To  be  consistent  In  our  viewpoint,  we  now  think  of  ^(y,t) 
as  the  density  of  neutrons  at  y  at  time  t.  Then  the  net 
neutron  current  density  l(y,t)  Is  provided  by  Fick's  Law,  in 
the  ordinary  diffusion  approach  [4o]  , 

(1)  i(y,t)  =  —  D  ^  j*(y,t). 

The  conservation  of  particles  (since  there  is  no  internal  pro¬ 
duction  when  k  »  1/2)  requires  in  any  interval  (a,b)  of  the 
rod 


(2) 


l(b,t)  -  i(a,t ) 


a 


Let  us  write,  for  the  net  current  emerging  from  our  rod  of 
length  x,  k(x,t).  Here,  again,  while  it  is  true  that 
k(x,t)  =  l(x,t)  we  choose  to  regard  the  x  in  the  function  k 
as  referring  to  the  length  of  the  rod.  Thus  k(x  +  A,t)  is  the 
net  current  emergent  from  a  rod  of  length  x  +  A,  source  2q(t) 
at  (x  +  A),  other  initial  and  boundary  conditions  being  as 
before . 

We  now  try  to  express  k(x  +  A,t)  in  terms  of  k(x,t). 
Applying  (2)  to  the  rod  of  length  x  +  A  we  find 

(3)  k(x  +  A ,  t )  -  l(x,t)  »  -  ^C/?X+A  /(y#t)dy, 

X 


or,  integrating  over  time, 

(4)  K(x  +  A,  t )  -  l(x,t)  =  “,/>)X+A  ^(y»t)dy. 

L  x 
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We  now  seek  expressions  for  l(x,t)  and  jrf.  To  find 
l(x,t)  we  note  that  we  have  thus  far  disregarded  the  part  of 
the  rod  from  0  to  x.  By  the  continuity  conditions  imposed 
by  diffusion  theory,  we  know  that  l(x,t)  is  merely  the 
current  out  of  x  due  to  the  source  ^(x,t)  imposed.  Let  us 
suppose  that  a  steady  source  of  unit  strength  produces  a 
current  out  of  the  rod  of  p(x,t).  Then  a  source  ^(x,t) 
will  produce  an  integrated  current 

(5)  I(x,t)  =  /;t  p(x,t  -  z)^(x,z)dz. 

(This  is  Just  Duhamel's  Principle  [35]). 

As  yet  we  have  not  used  (l).  From  it  we  find 

(6)  ^(x,t)  =  k(x  +  A ,  t )  +  2q(t)  +  o(  A  )  . 

Substituting  (5)  and  (6)  in  (4),  we  obtain 

(7)  K(x  +  A ,  t )  -c/)t  P(x,t  -  z)j^  k(x  +  A ,  z )  +  2q(  z  )l  dz 

0 

-  2Aq( t )  +  o( A ) . 

But,  bv  Duhamel's  Principle, 

(8)  /;t  p(x,t  -  z  )  2q(  z  )dz  =  K(x,t). 

0 

Thus 

ji +  2q(t)  -  p<x>t  -  2)“(x.z)dZ. 


(9) 
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This  agrees  with  (44.8)  upon  identifying  k  with  DhQ(q) 
and  p  with  hQ(l),  the  factor  1/2  occurring  because  p 
is  the  current  due  to  a  unit  source,  while  hQ  is  obtained 
from  a  source  of  strength  2. 

It  is  clear  that 

(10)  K(x,0)  =  0. 

To  find  K(0, t )  we  note  from  (6)  that 

^(0,t)  =  0  =  k(  A ,  t )  +  2q(  t )  +  o(  A  ) , 
so  that  for  Q( t )  >  0, 

(11)  k(o, t )  =  -  00  . 

Clearly,  in  case  q  =  1 ,  we  have 

(12)  H  +  2  =  p(x,t  -  z)p(x,z)dz, 

P(x,0 )  =  0,  P(0,t )  =  -  00  . 
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Part  IV 

> 

RANDOM  WALK  AND  MULTIPLE  SCATTERING 

46.  Random  Walk 

We  now  wish  to  apply  invariant  imbedding  techniques  to  the 
study  of  random  walk  processes.  Subsequently,  we  shall  consider 
more  general  processes  of  this  nature,  equivalent  to  multiple 
scattering  processes. 

Consider  the  finite  one-dimensional  lattice  consisting  of 
the  integer  values  between  a  and  b  along  the  real  line,  as 
indicated  below. 


, - 1 - f - 1 - 1 - 1 

a  a+1  a+2  k  b—1  b 

Fig.  19 

A  particle  Jumps  from  lattice-point  to  lattice-point  in 
accordance  with  the  following  law.  When  at  k,  there  is  a 
probability  q(k)  of  moving  one  unit  to  the  right  and  a 
probability  p(k)  =  1  —  q(k)  of  moving  one  unit  to  the  left. 

We  wish  to  determine  the  probability  that  the  particle, 
starting  at  k,  hits  the  barrier  at  a  before  it  hits  the 
barrier  at  b.  The  process  ends  as  soon  as  the  particle  lands 
at  a  or  b--hence,  the  name  "absorbing  barriers." 
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This  is  an  extension  of  the  classical  "gambler's  ruin" 
process  in  which  p(k)  •  q(k)  ■  1/2  for  a  <  k  <  b.  This 
particular  case  can  be  treated  in  a  very  elegant  fashion  by 
■•ana  of  Wald’ a  "fundamental  identity."  See  [j6]  for 
generalizations  the  case  of  dependent  steps.  Another 
extension  is  the  "game  of  survival,"  [18]  . 

The  classical  approach,  also  based  upon  the  use  of 
x'wourrence  relations,  proceeds  as  follows,  [j>7]  .  Let 

(1)  u(k)  -  the  probability  that  a  particle  starting  at 

k  lands  at  a  before  it  lands  at  b. 

Then,  considering  what  happens  as  a  result  of  the  first  step, 
for  a  <  k  <  b, 

(2)  u(k)  -  p(k)u(k  -  1)  +  q(k)u(k  +  l), 
with  the  two— point  boundary  conditions 

(3)  u(a)  -  1,  u(b)  -  0. 

We  thus  face  the  problem  of  solving  a  system  of  linear 
equations,  something  we  wish  to  avoid  if  possible.  In  what 
follows  we  shall  attack  these  problems  by  a  quite  different 
method. 

47.  Invariart  Imbedding  Approach 

The  observation  that  the  desired  probability,  u(k),  is 


a  function  of  ♦‘he  endpoints  a  and  b,  as  well  as  of  k, 
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kjya  our  approach.  Hence,  we  should  write  u(k,a,b)  to 
signify  this  dependence. 

Invariant  Imbedding,  as  we  have  repeatedly  stated, 
capitalizes  upon  this  dependence.  In  place  of  considering 
a  and  b  to  be  constants,  we  consider  them  to  be  parameters 
of  equal  importance  with  k,  which  means  that  in  place  of 
treating  an  individual  random  walk  process,  we  investigate 
simultaneously  an  entire  family  of  processes.  The  individual 
problem  is  then  analyzed  in  terms  of  its  relation  to 
contiguous  members  of  the  family.  In  this  way,  we  hope  to 
construct  a  bridge  between  particular  elements  of  the  family 
of  quite  simple  analytic  structure  and  the  processes  of 
interest  to  us. 

Let  us  then  keep  one  endpoint  fixed,  say  b,  and  regard 
u  as  a  function  of  a  and  k.  To  make  this  dependence 
explicit,  let  us  introduce  the  new  function 

(l)  f(a,k)  *  the  probability  that  a  particle  starting  at  k, 

a  ^  k  ^  b,  will  reach  a  before  reaching  b. 

Clearly  f(a,k)  -  u(k;a,b). 

To  obtain  the  required  relation  between  contiguous  ele¬ 
ments,  we  use  the  simple  geometric  fact  that  a  particle 
starting  at  k,  and  moving  one  unit  in  either  direction  at 
each  stage,  must  hit  a  +  1  before  it  can  hit  a.  This  is 
analogous  to  the  stratification  technique  we  used  in  the 
discussion  of  neutron  transport  and  multiplication.  Having 
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reached  a  +  1,  the  particle  must  then  reach  a  before  b, 
if  this  is  to  be  the  ease  for  the  original  process. 

Translating  these  remarks  into  algebraic  relations, 
using  the  elementary  rules  of  probability  theory,  we  have 

(2)  f(a,k)  -  f (a  4  l,k)f(a,a  +  1). 

Iterating  this  relation,  we  have 

(3)  f (a,k )  =  f(a,a  +  l)f(a  +  l,a  +  2)*--f(k  -  l,k). 

an  interesting  representation  for  solutions  of  a  Jacobi 
system  of  equations  of  the  type  appearing  in  (57.2). 

We  have  deliberately  stressed  the  word  "elementary”  in 
the  foregoing  discussion,  since  we  can  employ  similar  methods 
based  upon  other  concepts  of  probability  to  discuss  other 
kinds  of  equations.  Por  other  methods  of  treating  Jacobi 
equations,  see  [38]. 

48.  The  Function  f(a,a  +  1) 

To  answer  our  original  question,  that  of  determining  the 
value  of  f(a,k),  it  remains  to  evaluate  the  functions  of 
one  variable 

(1)  g(a)  -  f (a,a  +  l), 

defined  for  a  £  b  -  1.  Clearly  g(b  —  1)  -  0. 

Reverting  to  the  description  of  the  original  process,  we 


obtain  the  relation 
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(2)  f(a,a  +  1)  =  p(a  +  1 )  +  q(a  +  l)f(a,a  +  2). 

Using  (47.2),  we  have 

(3)  f(a,a  +  2)  *  f(a  +  l,a  +  2)f(a,a  +  1). 

Combining  these  two  expressions,  we  obtain  the  recurrence 
relation 


Since,  as  noted  above,  g(b  -  1)  =  0,  we  have  a  simple 
Inductive  determination  of  u(a)  for  a  <;  b  -  1. 

49.  An  Alternative  Derivation 

The  foregoing  result  can  be  derived  In  a  way  which  empha- 
sizer  Its  physical  significance  and  Its  connection  with 
previous  work  on  neutron  diffusion.  Toward  this  end,  let  us 
consider  the  following  scattering  problem.  A  rod  extending 
from  a  to  b  has  the  property  that  If  a  particle  Is  at  the 
position  k,  there  Is  probability  p(k)  that  It  will  be 
scattered  to  k  —  1,  and  probability  q(k)  *=  1  —  p(k)  that 
It  will  be  scattered  to  k  +  1.  A  particle  Is  placed  at  the 
end  a,  and  we  wish  to  determine  the  probability  that  It  will 
be  "back-scattered"  (reflected)  from  a,  over  all  time,  rather 
than  be  " forward- scattered"  (transmitted)  through  the  end  b. 

de  Imbed  this  process  within  the  class  of  processes  with 
"trigger'  particles  placed  at  the  end  AL  of  rods  extending 
from  jt  to  b,  with  A  =  b,b  -  1,...,  and  then  write  a 
functional  equation  interconnecting  these  processes. 
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Let  us  define  the  function  g(a)  directly 

(1)  g(a)  c  the  probability  that  over  all  time  a  particle 

at  a  +  1  will  be  backscattered  to  a  by  the 
rod  extending  from  a  +  1  to  b,  rather  than 
be  forwardscattered  at  b  from  the  rod. 

Observe  next  that  with  the  particle  initially  at  the  position 
a  +  1,  there  Is  probability  p(a  +  1)  that  it  will  be 
scattered  directly  to  the  point  a.  On  the  other  hand,  if  it 
is  scattered  initially  to  the  right  to  (a  +  2),  then  by 
definition  there  is  probability  g(a  +  l)  that  it  will 
eventually  be  back-scattered  from  the  rod  (a  +  2,b)  to  the 
point  a  ¥  1,  from  which  it  may  be  scattered  to  the  point  a. 
Should,  however,  it  once  again  be  scattered  to  the  right,  to 
the  point  (a  +  2),  there  is  once  again  probability  g(a  +  1) 
that  it  will  eventually  reach  the  point  a  +  1,  and  so  on. 

In  this  way  we  see  that 

(2)  g(a )  -  p(a  +  1)  +  q(a  +  l)g(a  +  l)p(a  +  1) 

+  q(a  +  1 )g( a  +  l)q(a  +  l)g(a  +  l)p(a  +  l)  +  •••, 

which,  upon  summing  the  geometric  series  on  the  right  hand 
side,  leads  to  the  equation 


(3) 


*<»> 


17* 
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This  Id  Equation  (48.4). 

The  method  that  we  have  used  In  deriving  Equation  (48.2) 
of  (3)  corresponds  abstractly  to  the  method  used  by 
Ambarzumlan  [29]  in  discussing  diffuse  reflection  from  a  foggy 
medium,  a  process  we  shall  discuss  below,  and  to  the  method 
used  above  in  handling  some  neutron  transport  processes. 

Similarly,  the  discussion  in  ^47  may  be  reinterpreted  to 
yield  the  flux  scattered  from  the  end  a  of  a  rod  as  a  result 
of  an  internal  source  of  particles. 

50.  Expected  Sojourn 

Let  us  now  introduce  the  function 

(1)  w(a,k)  =  the  conditional  expected  time  required  for 

the  particle  to  reach  a  before  reaching 
o,  starting  from  k,  assuming  a  unit  step 
takes  unit  time. 

By  this  we  mean  the  expected  time  required  to  reach  a,  under 
the  assumption  that  a  is  reached  before  b.  Then,  the  same 
reasoning  as  above  yields  the  relation 

(2)  w(a,k)  =  w(a  l,k)  +  w(a,a  +  1). 

To  obtain  an  analytic  expression  for  w(a,a  +  1 ) ,  we  combine 
the  two  expressions 

(3)  w ( a , a  +  1)  =  p(a  +  1)  f  q(a  +  l)[w(a,a  +  2)  +  l]  , 


w(a,a  +  2)  «=  w(a  +  l,a  +  2)  +  w(a,a  +  1). 
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The  result  Is 


(4) 


"(a>a  +  -  FTTTT}  +  p[a  7 4 j  “(a  +  x'a  +  2)- 


Iteration  yields  the  infinite  series 


(5) 


w(a’a  +  x)  "  FtTTTT  +  Vll 


_ q (a  +  1 )q(a  +  2 ) _ 

+  p  ( a  +  3  )  P  ( a  +  2 )  p  ( a  '+  1 7 


with  the  convention  that  the  series  terminates  If  b  is  finite. 


51.  Characteristic  Functions 

As  we  have  pointed  out  above,  whenever  an  expected  value 
can  oe  determined,  the  same  techniques  yield  relations  for 
generating  functions.  Let 

(1)  y(a,k,s)  -  a(elBZ), 

where  z  =  z(a,k)  is  the  random  variable  equal  to  the  time 
spent  by  the  particle  in  going  from  k  to  a,  without  ever 
hitting  b. 

Since 


(2)  z(a,k)  =  z(a  +  1 , k )  +  z(a,a  +  1), 
we  have  the  functional  equation 

(3)  y(a,k,s)  =  y(a  +  l,k,s)y(a,a  +  l,s). 


It  Is  now  not  hard  to  show  that 
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(4)  y(a,a  4  1)  ■»  p(a  4  l)els  4  q(a  4  l)elsy(a,a  4  2), 

y(a,a  4  2)  *  y(a  4  l,a  4  2)y(a,a  4-  l). 


From  these  we  obtain  the  recurrence  relation 

is 


(5) 


y(a,a  4  l)  - 


P  (?  +  !>; 

is 


1  -  q(a  4  l)e  y(a  4  l,a  4  2) 


which  can  be  used  to  obtain  higher  momenta,  ] . 


52.  More  Qeneral  Random  Walk  Processes 

Similar  methods  can  be  applied  to  random  walks  which 
allow  steps  of  different  units  at  each  stage.  Abstractly, 
these  methods  will  be  equivalent  when  vectorwnatrix  notation 
1 s  introduced. 


53,  Multiple  Scattering 

In  this  part  of  the  paper,  we  wish  to  consider  two 
processes  of  particular  interest.  The  first  concerns  a  one- 
dimensional  random  walk  in  which  the  energy  of  the  particle 
changes  as  a  result  of  its  wandering,  while  the  second  per¬ 
tains  to  a  two-or-three-dimensional  random  walk  process  in  which 
the  direction  of  motion  changes  at  a  result  of  each  collision. 

Both  of  these  can,  in  discrete  form,  be  considered  to  be 
particular  cases  of  a  process  of  the  following  type: 

■  "A  particle  in  state  i,  1  ■  1,2,...,N,  can  occupy  any 
of  the  lattice  points  k  between  a  and  b.  Let 
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(1)  p^j(k)  =  Pro'oabllity  that  a  particle  at  k  in 

state  i  will  go  one  unit  to  the  left,  and 
arrive  in  state  j,  i , J  =  1,2, ...,N, 

q^j(k)  =  the  probability  that  a  particle  at  k  in 
state  i  will  go  one  unit  to  the  right, 

« 

and  arrive  in  state  J,  1,J  =  1 , 2, . . . ,N . " 

2 

Let  us  then  define  the  N  functions 

(2)  u^j(a,k)  =  the  probability  that  a  particle  starting 

at  k  in  state  i  will  hit  a  in  states 
J  before  reaching  b  in  any  state. 

Proceeding  as  in  the  foregoing  sections,  we  obtain  the 
relation 

(3)  ui j(a»k)  Z  ulm(a  a  l,k)umJ(a,a  a  l). 

Consequently,  if  we  introduce  the  matrix  function 

(4)  U(a,k)  -  (ui ^(a,k) ) , 
we  derive  the  basic  relation 

(5)  U(a,k )  =  U( a  +  l,k)u(a,a  +  1), 

the  analogue  of  (47.2).  Once  again,  we  see  that  the  problem 
has  reduced  to  a  determination  of  a  function  of  a  alone, 

U(a,a  +  1 ) . 

_ 

It  is  convenient  here  to  reverse  the  index  order  as  com¬ 
pared  to  previous  usage  (see  §8). 
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54 .  Determination  of  U(a,a  +  1) 

Proceeding  as  before,  we  have 

(1)  ulj(a'a  +  x)  =  pij(a  +  x)  +  z  Qlrn(a  +  1 ^UraJ ^ a,a  +  2)» 

or.  If  0(a)  -  (Ujj (a,a  +  1))#  P(a)  =  (p1J(a)),  Q( a )  =  (q^a)), 

(2)  G(a)  -  P(a  +  1)  +  Q(a  +  l)U(a,a  +  2). 

Using  (53.2)#  we  obtain  the  relation 

(3)  U(a,a  4.  2)  =  U( a  +  l,a  +  2)U(a,a  4-  1). 

Using  this  In  (2),  we  have 

(4)  0(a)  =  P(a  +  1)  +  Q( a  +  l)0(a  +  l)G(a), 
or 

(5)  G( a )  =  [i  -  Q(a  +  1  )Q( a  +  l)]_ip(a  +  1). 

Since  G(b  —  1)  =  0,  once  again  we  have  a  direct  Iterative 
technique  for  determining  G(a),  and  thus  U(a,x). 

55.  Discussion 

It  Is  clear  that  In  the  same  way  we  can  obtain  multidim¬ 
ensional  analogues  of  the  results  for  the  expected  sojourn  and 
generating  function. 

Turning  from  the  analytic  aspects,  let  us  examine  the 
computational  aspects.  Approaching  the  problem  along  conven¬ 
tional  lines,  we  obtain  a  system  of  linear  equations  of  the  form 
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(1)  U  (k)  *  z  Plmu  (k  -  1)  +  q  (k)u  (k  +  i). 

0  m  °  m  0 

Keeping  J  fixed,  we  have  a  system  of  order  N(b  —  a).  If 
N  *  10,  and  b  —  a  =  100,  this  is  order  1000,  a  respectable 
and  even  formidable  number,  even  in  the  light  of  modem  devices 
Using  the  technique  described  above,  the  solution  is  made 
to  depend  upon  the  inversion  and  repeated  multiplication  of 
10  x  10  matrices,  a  complicated,  but  far  more  feasible  process 

56.  Time-dependent  Processes 

Let  us  now  consider  the  one-dimensional  random  walk  pro¬ 
cess  in  which  each  stage  consumes  one  time  unit.  Let 

(1)  u(a,k,t)  =  the  probability  of  going  from  k  to  a 

in  time  t,  without  ever  hitting  b. 

As  we  shall  see,  we  derive  equations  completely  analogous  to 
those  exhibited  above  for  the  generating  function 

00  t 

(2)  F(a,k,t,r)  s  F(a,k)  =  Z  u(a,k,t)r\ 

t  =0 


As  in  the  previous  sections,  we  obtain  the  fundamental 
relation 

t 

(3)  u(a,k,t)  =  Z  u(a  +  l,k,s)u(a,a  +  l,t  -  s), 

s=0 


and  once  again  derive  the  fundamental  relation 

(4)  ‘  F(a,k)  =  F(a  +  l,k)F(a,a  +  l). 
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Furthermore, 

(5)  u(a,a  +  l,t)  =  p(a  +  l)b(l,t)  +  q(a  +  l)u(a,a  +  2,t  —  1), 
for  t  ^  1,  where 

(6)  b(l, t )  =1,  t  -  l, 

-  0,  t  /  1. 

Thus,  multiplying  (5)  by  rfc  and  summing  over  t, 

(7)  F(a,a  +  1)  =  p(a  +  l)r  +  q(a  +  l)rF(a,a  +  2). 


From  here  on  the  argument  proceeds  as  before.  The  final 
result  Is 


F(a,a  +  1)  =  y-^r 


q(a  + 


a  +  1  )r 

‘l )rF(a  +  l,a  +  2 


Equation  (7)  can  be  derived  directly  making  use  of  the 
properties  of  the  process  and  of  generating  functions.  It 
should  be  compared  with  (51*5) • 
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Part  V 

RADIATIVE  TRANSFER 


57.  Introduction 

In  these  last  few  sections  we  shall  discuss  a  problem 
arising  in  the  field  of  radiative  transfer.  Abstractly,  such 
cases  are  equivalent  to  appropriate  neutron  transport  problems, 
and  therefore  the  material  appearing  here  could  very  well  have 
been  placed  in  earlier  sections.  However,  since  it  was  in  the 
solution  of  problems  of  this  genre  that  Arabarzumlan  [29]  first 
successfully  used  his  invariance  principle  and  here,  too,  that 
Chandrasekhar  developed  his  extensive  genbrallzation  30J  it 
seems  fitting  that  we  leave  these  problems  In  their  original 
setting.  Our  discussion  follows  that  of  [42]  in  which  the 
principle  of  invariant  imbedding  was  first  sketched. 

58.  The  Physical  Model 

Assume  that  parallel  rays  of  light  of  uniform  Intensity 
are  incident  on  an  inhomogeneous  slab  composed  of  a  substance 
which  absorbs  and  scattero  light.  Our  objective  is  to  deter¬ 
mine  the  Intensity  of  the  diffusely  reflected  light  as  a 
function  of  the  Incident  light,  the  properties  of  the  slab, 
and  the  angle  of  the  emerging  rays  . 
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0 

Pig.  20 

Cross  section  of  slab  with 
Incident  and  reflected  rays. 

We  shall  assume  the  slab  has  the  following  absorption  and 
scattering  properties. 

(1)  In  traversing  a  distance  d  In  the  slab  a  portion 
of  a  beam  I  Is  reduced  to  Intensity  1(1  -  ad)  +  o(d).  A 
fraction  7\  of  the  Intercepted  beam  is  reradiated,  while  a 
fraction  1  -  7\  Is  permanently  lost  (absorbed).  The 
quantities  "A  and  a  may  be  dependent  upon  the  distance  from 
an  edge  of  the  slab. 

(2)  Radiation  is  scattered  isotropically. 

In  the  light  of  our  previous  work,  it  is  most  convenient 
to  consider  the  light  radiation  as  composed  of  photons; 
particles,  then,  which  behave,  for  our  purposes,  Just  as 
neutrons  do. 

We  must  note,  too,  that  In  reality  not  only  the  angles 
Y  and  0  (Pig.  20)  arise  but  also  the  corresponding 
azimuthal  angles.  The  latter  may  be  neglected  in  our  analysis 
because  of  the  symmetry  of  the  problem. 
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We  define  a  reflection  coefficient  R(x,vf,0),  giving  the 
intensity  of  radiation  reflected  in  direction  0  per  unit  area 
on  the  face  of  the  slab  due  to  a  beam  of  unit  intensity  inci¬ 
dent  at  angle  'f,  the  thickness  of  the  slab  being  x. 

Otherwise  put,  R  is  the  number  of  photons  per  unit  time  out 
of  a  unit  area  on  the  slab  at  angle  0,  due  to  a  unit  flux 
(one  photon  per  unit  area  per  unit  time)  impinging  at  angle  Y* 
It  is  now  easy  to  write  down  the  equations  for  R,  either 
using  directly  the  fundamental  principles  we  have  developed,  or 
applying  the  general  flux  equation  (26.2).  We  find,  recalling 
that  integrations  over  the  azimuthal  angles  are  necessary, 
even  though  R  is  independent  of  them, 


(1) 


2R 

3x 


coTU  + 

R(x,f • ,0)sin  f 'df ' 


+  a(x)7\(x)  n 
2  °X) 


+  a(x)7\(x)v  P 
Vq 


t/2 


R(x,4\e') 


sin  9 1 
cos  0' 


d0‘ 


cos  6'  .etalnf'd'p. 


59.  Comparison  with  the  Results  of  Ambarzumlan 

Equation  (58.1)  is  considerably  more  general  than  the 
original  result  of  Ambarzumlan,  since  that  writer  considered 
a  semi— inf inlte  slab,  with  a  »  1  and  1\  constant.  Our 
result  should  hence  reduce  to  his  when  we  eliminate  the  x 
dependence  and  set  3 R/dx.  -  0.  The  reader  will  find,  however, 
that  the  equations  still  differ  considerably. 
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The  reason  for  this  apparent  discrepancy  lies  in  the  way 
we  have  chosen  to  measure  flux  throughout  this  paper.  While 
classically  flux  is  measured  in  terms  of  the  number  of  particles 
per  unit  time  crossing  a  unit  area  normal  to  the  direction  of 
the  particle,  we  have  chosen  instead  to  talk  in  terms  of  the 
number  per  unit  time  through  a  unit  area  on  the  (geometric) 
surface  through  whxch  the  particles  are  passing. 

A  bit  of  philosophy  may  be  appropriate  .  The  transport 
equation  for  the  flux  of  particles  internal  to  a  body 
ordinarily  is  quite  independent  of  the  boundaries  of  the 
medium  Itself.  The  geometry  is  brought  in  through  the 
auxiliary  boundary  conditions.  Thus  a  person  in  outer  space, 
with  no  frame  of  reference,  would  rather  naturally  measure 
flux  in  the  classical  way. 

Tne  equations  of  invariant  imbedding,  however,  depend 
deeply  and  Inherently  upon  the  boundaries  of  the  body  under 
consideration.  It  therefore  seems  natural  to  define  the  flux 
with  direct  reference  to  those  boundaries  .  We  have  found 
this  the  easier  way  conceptually. 

In  any  event  it  is  possible  to  convert  from  one  definition 
of  flux  to  the  other  without  great  effort.  Consider  the  case 
of  §58,  and  let  us  define  l(x,"f,0)  as  the  reflected  intensity 
in  the  classical  senBe:  l(x,f,0)  is  the  number  of  photonB 
reflected  from  x  travelling  in  direction  0  per  unit  time 
through  a  unit  area  normal  to  that  direction  due  to  one  photon 
Incident  on  x  per  unit  time  per  unit  area  normal  tc  the 
direction  vf.  Then  it  is  easy  to  see  that 
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(i)  i(x,t,e)  »  R(*/f,6)(££L^). 


Using  the  foregoing  transformation,  we  obtain  Ambarzumlan ' b 
equation . 
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? art  VI 
SUMMARY 

60.  Review  of  Basic  Techniques 

Having  covered  some  quite  diverse  parts  of  mathematical 
physics,  we  feel  that  It  is  important  to  state  to  the  reader 
what  our  basic  ideas  and  objectives  have  been,  and  what  have 
been  the  methods  we  have  employed  toward  obtaining  thase  goals. 

We  start  with  the  fact  that  any  physical  process  can  be 
described  in  a  variety  of  different  ways,  leading  to  a  number 
of  different  analytic  paraphrases.  As  soon  as  thl3  mo3t 
important  fact  1s  accepted,  then  necessarily  the  premise  must 
be  accepted  that  some  descriptions  will  be  significantly 
better  than  others  for  the  3tudy  of  particular  properties  of 
a  process . 

We  have  hinged  our  description  of  physical  processes  upon 
the  invariance  concept .  By  this  we  mean  that  we  have  con¬ 
sistently  introduced  state  variables  and  written  our  equations 
in  such  a  way  as  to  stress  the  idea  that  any  individual  pro¬ 
cess  is  to  be  considered  as  a  member  of  a  family  of  related 
processes.  The  advantage  to  be  gained  from  this  point  of  view 
resides  in  the  common  observation  that  the  properties  of  a 
particular  member  of  a  set  can  often  be  easily  understood  in 
terms  of  the  properties  of  contiguous  elements,  although  often 
quite  puzzling  to  comprehend  in  isolation.  This  principle  of 
continuity,  one  of  the  most  powerful  and  versatile  tools  in 
the  mathematicians  hope  chest,  is  basic  also  in  the  biological 


world . 
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The  usual  equations  of  mathematical  physics  arise  from 
the  application  of  imbedding  techniques,  by  means  of  the 
Introduction  of  fluxes  at  arbitrary  points,  by  means  of  the 
introduction  of  time,  and  so  on .  By  introducing  other  para¬ 
meters  of  significance  and  applying  invariance  principles  in 
a  different  fashion,  we  have  obtained  new  equations  which 
possess  certain  computational  and  analytic  advantages  over 
the  classical  formulations.  We  have  indicated  how  one  can 
pass  back  and  forth  from  one  type  of  equation  to  the  other. 

One  of  the  major  difficulties  of  the  classical  approach 
lies  in  the  fact  that  it  leads  to  boundary  value  problems 
which  in  turn  lead  to  Fredholm  integral  equations.  Ultimately 
for  the  solution  of  these  problems  in  numerical  terms  we  are 
forced  to  the  solution  of  large  systems  of  linear  equations. 
TYiis  is  a  most  subtle  and  difficult  problem  and  one  for  which 
modern  digital  computers  are  not  well  suited. 

The  new  approach  presented  here  leads  to  the  solution  of 
initial  value  problems,  not  necessarily  in  time  but  in  other 
meaningful  physical  parameters,  and  computationally  to  the 
iteration  of  nonlinear  transformations.  This  latter  is  a 
task  well  designed  for  the  digital  computer.  Finally,  we 
note  that  our  formulation  very  often  corresponds  to  the  way 
the  data  is  obtained  experimentally. 

Let  us  note  in  passing  that  there  are  available  other 
techniques  for  converting  boundary  value  problems  involving 
Fredholm  type  Integral  equations  to  Initial  value  problems 
relying  upon  Volterra  type  integral  equations.  Generalizing 
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results  of  Holmgren  and  Levi,  Muntz  discuBsed  how  this  could 
be  done  for  the  heat  equation,  while  Milgram  and  Roaenbloom 
used  a  different  device  for  treating  generalized  potential 
theory#  involving  the  harmonic  integrals  of  Hodge.  Our 
methods  are  quite  different  from  theirs,  since  ours  always 
involve  a  variation  of  the  domain,  while  theirs  keep  the 
domain  fixed . 

Secondly,  let  us  point  out  that  the  Riccati  equations 
which  are  characteristic  of  our  approach  can  be  obtained  in  a 
number  of  different  ways  in  dealing  with  ordinary  second  order 
differential  equations.  The  derivation  we  employ  is  quite 
different  from  that  obtained  by  a  simple  change  of  a  dependent 
variable,  and  different  also  from  that  resulting  from  dynamic 
programming.  Our  methods  lead  naturally  to  the  generalized 
Riccati  equations  corresponding  to  partial  differential 
operators  . 

The  functions  with  which  we  deal,  the  reflection  and 
transmission  functions,  appear  to  be  basic  functions  of  analy¬ 
sis.  IMs  is  to  be  expected  since  they  represent  fundamental 
physical  quantities.  Using  these  functions  we  obtain  a  new 
approach  .to  the  problems  of  existence  and  uniqueness  of 
solutions  of  the  classical  equations,  and  new  representations 
for  the  solutions  of  these  equations.  What  is  most  important 
is  that  these  methods  are  independent  of  characteristic  value 
and  spectral  theory.  Results  of  this  type  will  appear  in  the 


near  future . 
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Although  we  have  from  time  to  time  mentioned  various 
computational  advantages  of  our  procedures,  we  have  not 
Included  any  calculations  in  this  pape1^  for  several  reasons. 
In  the  first  place  the  paper  has  already  assumed  a  certain 
unwieldy  length.  Par  more  Important  Is  the  fact  that 
numerical  solution  of  significant  problems  introduce  a  number 
of  nontrivial  questions,  regardless  of  the  luethoda  that  s 
used.  Any  who  have  engaged  In  the  computational  solution  of 
equations  will  sadly  testify  to  this.  Consequent i„  ,  we  feel 
that  it  la  better  to  leave  the  calculation*?  for  a  separate 
study,  devoted  merely  to  thi8  aspect. 
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